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ABSTRACT 

We present three-dimensional (3D) simulations of supernova explosions of nonrotating stars, triggered by the delayed neutrino- 
heating mechanism with a suitable choice of the core-neutrino luminosity. Our results show that asymmetric mass ejection caused by 
hydrodynamic instabilities can accelerate the neutron star (NS) up to recoil velocities of more than 700kms~^ by the "gravitational 
tug-boat mechanism", which is sufficient to explain most observed pulsar space velocities. The associated NS spin periods for our 
nonrotating progenitors are about 100 ms to 8000 ms without any obvious correlation between spin and kick magnitudes or directions. 
This suggests that faster spins and a possible spin-kick alignment might require angular momentum in the progenitor core prior to 
collapse. Our simulations for the first time demonstrate a clear correlation between the size of the NS kick and anisotropic production 
and distribution of heavy elements created by explosive burning behind the shock. In the case of large pulsar kicks the explosion is 
significantly stronger opposite to the kick vector. Therefore the bulk of the explosively fused iron-group elements, in particular nickel, 
is ejected mostly in large clumps against the kick direction. This contrasts with the case of low recoil velocity, where the nickel-rich 
lumps are more isotropically distributed. Explosively produced intermediate-mass nuclei heavier than ^^Si (like "^^Ca and "^"^Ti) also 
exhibit a significant enhancement in the hemisphere opposite to the direction of fast NS motion, while the distribution of ^^C, ^^O, and 
^^Ne is not aff'ected, and that of ^^Mg only marginally. Mapping the spatial distribution of the heavy elements in supernova remnants 
with identified pulsar motion may off'er an important diagnostic test of the kick mechanism. Diff'erent from kick scenarios based on 
anisotropic neutrino emission, our hydrodynamical acceleration model predicts enhanced ejection of iron-group elements and of their 
nuclear precursors in the direction opposite to the NS recoil. 

Key words. Supernovae: general — pulsars: general — stars: neutron — hydrodynamics — neutrinos 
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Young pulsars are measured to possess high sp ace velocities 
with an average value of several lOOkms"^ (e.g., Arzoumanian 



let al.|2QQ2||Hobbs et al.|2Q05l|Faucher-Giguere & Kaspi|2Q06| ) 
Some of them move with an estimated transverse velocity even 
higher than 500kms"\ like e.g., the compact stellar remnant 
RX J0822-43Q0 in Pu ppis A with a velocity of 672 ± 1 15 km s"^ 
( [Becker et al.||2Q 12 |), and PSR B 1508+ 55 with a velocity of 



083^^^^ km s 



(Chatterjee et al. 



2005 1). Such velocities are 
he disruption of close 



much are too high to be explained by 
binary systems in the supernova (SN) event that gave birth to 
the neutron star (NS). Further evidence for natal kicks and as- 
sociated NS spin-up can be deduced from orbital and spin pa- 
rameters and kinematic informatio n of doubl e NS systems (e.g., 
^ong et a l. 2010; Farr et al. 20lTt[LapQQT] ). 



Although only one percent anisotropy of the several 10^^ erg 
of gravitational binding energy released through neutrinos in the 
NS formation event would be sufficient to account for a NS re- 
coil velocity of about 300kms"\ it turns out that an asymme- 
try of this size or larger is extremely difficult to obtain. Super- 
strong (^ 10^^ G) internal, ordered magnetic fields, i.e., fields 
with a very strong dipolar component, would be needed in the 
NS. Such fields could impose a direction dependence on the 
matter-interactions of active neutrinos in the de nse NS core (e.g., 
Bisnovatyi- Kogan 1996; Lai & Qian 1998 ; Sagert & Schaff^ner-] 
Bielich 2008|), or they :ould create an asymmetry of resonant 
Mikheev-Smimov-Wolfenstein flavor transformations between 
active neutrinos and sterile neutrinos with masses of several 
keV (e.g., [Fryer & Kusenko"'2006^, 'Kusenko "2009^, 'Kishimotoj 
2011). Also a radiatively driven magnetoacoustic instability, 
termed neutrino-bubble instability, has been proposed as a pos- 



sible cause of sizable neutrino emission anisotropics ( Socrates 



Numerous models for producing the kicks have been pro- 
posed, either attempting to explain the NS recoil by anisotropic 
neutrino emission during the Kelvin-Helmholtz cooling of the 
nascent remnant or by asymmetric mass ejection in the SN ex- 
plosion (cf. [Lai et al.| 2001 , for a review). Both scenarios con- 
stitute fundamental possibilities that exploit momentum conser- 
vation in the transition from the progenitor star on the one hand 
to the relic compact object, expelled gas, and -0.2 M© of mass- 
equivalent neutrino loss on the other hand. The challenge, how- 
ever, is connected to the explanation of the exact physical reason 
of the necessary asymmetry. 



et al.||2005| ). However, all these scenarios require the favorable 
combination of several unestablished ingredients. 

In contrast, NS kick theories based on anisotropic mass ejec- 
tion during the explosion can be motivated by the observational 
fact that core-collapse SNe and their gaseous remnants exhibit 
large-scale asphericity essentially in all observed cases. While 
the fundamental requirement of the "hydrodynamical" recoil 
mechanism is thus well consolidated, the exact origin and size 
of the mass ejection asymmetry needs to be worked out and re- 



quires hydrodynamical modeling for quantitative answers. Janka 



& MUller] ( |1994| ), using numerical simulations and analytic es- 
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timates, showed that the interaction of the proto-neutron star 
(PNS) with the violent mass flows in the neutrino-heated post- 
shock layer prior to and around the onset of the SN explosion 
cannot explain NS velocities of many 100 km s"^ This conclu- 
sion is supported by more recent three-dimensional (3D) simula- 
tions by |Fryer & Young ( 2007| ). The mentioned works relied on 
asymmetries that developed only after core bounce by neutrino- 
driven convection and the standing- accretion shock instability 
(SASI; [Blondin et aL][2QQ3| |Foglizzo||2Q02l ). These hydrody- 
namic instabilities can grow from random seed perturbations, 
which can be expected to be present with small amplitudes in any 
convectively stirred stellar environment. In contrast, Burrows & 
Hayes| (p"996 ) assumed the existence of a large, dipolar asymme- 
try already in the pre-collapse stellar core as the relic of convec- 
tive shell burning during the late evolution stages of the progen- 
itor star. With the density accordingly reduced in a wide polar 
funnel, they obtained a pronounced pole-to-pole asymmetry of 
the explosion, because the ejecta could accelerate more readily 
to higher velocities on the side of the low-density funnel. As a 
consequence, the NS received a kick of more than 500 km s"^ 
opposite to the direction of the stronger explosion. Although 
the assumed density reduction might be motivated by the obser- 
vation of vigorous dynamical interaction between difl'erent burn- 
ing shells in combination with large (low-order mode) deviations 
from spherical symmetry in two-dimensional (2D) stellar calcu- 
lations for a stage close to core collapse ( A rnett & Meakin|201 
the amplitude and structure of such density and velocity per- 
turbations in pre-collapse stellar cores are still highly uncertain. 
The determination of these perturbations requires full-sphere 3D 
simulations over much longer periods of pre-collapse evolution 
than the existing 2D models. 

IScheck et"aL] ( |2004| |2006| ) therefore also started with 
minimal and least predetermining assumptions about inhomo- 
geneities in the progenitor before collapse. Like| Janka & MuUer] 
( 1994| ) they applied small random perturbations (typically 0. 1 %- 
1% of the density and/or velocity) to trigger the development of 
large-scale, nonradial flows in the postshock accretion layer, but 
in contrast to the previous work they could follow the asymmet- 
ric, neutrino-powered SN blast wave well beyond the onset of 
the explosion for evolution times of at least one second (corre- 
sponding to shock radii of -10,000 km and more). Their results 
for a large set of 2D simulations demonstrated that due to the 
action of neutrino-driven convection and low-mode SASI activ- 
ity, the postshock shell can be expelled with sufficiently large 
asymmetry to accelerate the NS over a timescale of seconds to 
velocities up to ~ 1000 km s"^ and more. The NS recoil devel- 
ops in the direction opposite to the side of the strongest ex- 
plosion. It is mediated mainly by the gravitational forces be- 
tween the compact remnant and the most slowly moving, mas- 
sive ejecta "clumps", whose radial propagation lags behind the 
rest of the expelled postshock material. The long-distance cou- 
pling by gravitational forces can still act on the NS long after 
hydrodynamical interactions have ceased. The influence of long- 
range gravity is therefore much more efficient in accelerating the 
NS into a well defined direction than the random kicks exerted 
by convective downdrafts during the accretion phase. The im- 
pacts of accretion flows bufl'et the accretor in varying directions. 
Such a random- walk like process is inefl'ective in producing high 
recoil velocities. Recent 2D results of ^Nordhaus et al.| ( |201Q I 
2012| lend support to the " gravitational tug-boat mechanism']^ 
and' Wong wathan arat et al.|p010b| ) presented a first, small set of 
3D simulations which showed that NS recoil velocities of at least 



-500 km s ^ can also be obtained in the absence of the artificial 
contraint to axisymmetry associated with 2D modeling. 

In this paper we discuss results of a significantly extended set 
of 3D simulations for nonrotating 15 M© and 20 M© progenitor 
stars, some of which were continued for nearly 3.5 seconds be- 
yond bounce to follow the saturation of the NS kick velocity. As 
in the previous works by |Scheck et al.| ( 2004, 2006) and j Wong- 1 
[wathanarat et al.|p010b| ), we initiate the neutrino-driven explo- 
sions artificially by imposing suitably chosen neutrino luminosi- 
ties at the inner grid boundary. The latter is placed well inside 
the neutrinosphere (at neutrino optically depths between ~10 and 
some 100), and neutrino eff'ects on the computational grid are 
treated by an approximate description of the neutrino transport. 
We also present estimates of the spins that are acquired by the 
NS. Moreover, we argue and demonstrate by results that the ex- 
plosive production of heavy elements during the first seconds of 
the SN blast becomes highly anisotropic in the case of large NS 
kicks. In particular nuclei between silicon and the iron group are 
ejected in considerably larger amounts in the hemisphere point- 
ing away from the NS velocity. We therefore predict that the 
bulk of such nucleosynthetic products should possess a momen- 
tum opposite to that of a high- velocity NS. Confirming such a 
hemispheric momentum asymmetry in SN remnants with deter- 
mined NS velocities would provide extremely valuable informa- 
tion about the underlying NS kick mechanism. In particular, it 
might allow one to discriminate our gravitational tug-boat mech- 
anism from scenarios like the neutrino-driven kicks discussed by 
Fryer & Kusenko| ( [2006| ), who expect the stronger explosion (and 



thus enhanced explosive nucleosynthesis) in the direction of the 
NS motion. High-resolution spectral investigations of young SN 
remnants like Pupp is A (Katsuda et al. 2008, 2010), Gil. 2 -0.3 
([Moon et al.||2009p, and Cassiopeia A (isensee et al.|2010t|De 



[Laney et a l.|2010[ |Rest et al.|2011t [Hwang & Laming|2012| ) in 

different wavebands may offer a promising perspective. 

The outline of our paper is as follows. In Sect. |2] we de- 
scribe our numerical setup, employed methods, and the com- 
puted models. In Sect. [3] we present our results for NS kicks, 
discuss their connection to blast- wave asymmetries (Sect. |3.2| ), 
give a detailed description of the physics of the NS acceleration 
mechanism (Sect. |3.3| ), provide analytic estimates of the achiev- 
able NS velocities (Sect. |3.4| ), and address possible progenitor 
dependences (Sect. |3.5| ) and the contribution from anisotropic 
neutrino emission (Sect. [T6| ). In Sect. |4] we evaluate our mod- 
els concerning the spin-up of the nascent NSs, and in Sect. [5] we 
discuss our expectations for the connection between NS kicks 
and observable asymmetries of the explosive nucleosynthesis of 
heavy elements. Finally we summarize and draw conclusions in 
Sect. [61 



2. Numerical methods and setup 

2.1. Code, grid, and boundary conditions 

Our numerical code is based on the explicit finite- volume, Eule- 
rian, multi-fluid hydrodynamics code Prometheus (Fryx ell et al.| 
|1991[[MUller et al.fl991 a b). It integrates the multidimensional 
hydrodynami c equations using dimensional splitting fo llowing 
Strangl('1968'), piecewise parabolic reconstruction (PPM; Colella 



We adopt this term following a suggestion by Jeremiah Murphy. decoupling ( |Quirk| 



& Woodw ard 1984 ), and a Riemann solver for real gases ( [Colella 
& Glaz|T98 5 ). Inside grid cells with strong grid- aligned shocks, 
fluxes computed from the Riemann solver are replaced by the 
AUSM-F fluxes of | Liou ( |1996| ) in order to prevent odd-even 
19941). The code treats advection of nu- 
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Fig. 1. Axis-free Yin- Yang grid plotted on a spherical surface demon- 
strating symmetric geometry of the Yin (red) and Yang (blue) regions. 
Both grids are the low-latitude part of a spherical polar grid, and there- 
fore contain grid cells which are almost equidistant in angular direc- 
tions. 



clear species by e mploying the Consistent Multi-fluid Advection 
(CMA) scheme of Plewa & Muller|( [T999 ^. 

The code employs an axis-free overlapping "Yin- Yang" grid 
( Kageyama & Sato|2004| ) in spherical polar coordinates, which 
was recently implemented into Prometheus, for spatial dis- 
cretization (^^^^^^^^^TeF^jSoTO^. Figure [T] shows the 
Yin- Yang grid plotted on a spherical surface. The Yin- Yang grid 
helps relaxing the restrictive CFL-timestep condition and avoid- 
ing numerical artifacts near the polar axis. Our standard grid 
configuration consists of 400(r) x 41(6) x 137(0) x 2 grid ceUs, 
corresponding to an angular resolution of 2° and covering the 
full 471 solid angle. The grid configuration with reduced angular 
resolution consists of 400(r) x 20(0) x 56(0) x 2 zones, corre- 
sponding to an angular resolution of 5°. The radial grid has a 
constant spacing of 0.3 km from the inner grid boundary 7?ib up 
to a radius of 80 km for the progenitors W15, N20, and B 15, and 
115 km for the progenitor L15 (see Sect. |2.5| ). Beyond this radius 
the radial grid is logarithmically spaced. The outer grid bound- 
ary 7?ob is at 18000-20000 km, which is sufficient to prevent the 
supernova (SN) shock from leaving the computational domain 
during the simulated time. Hydrostatic equilibrium is assumed 
at the inner grid boundary 7?ib, while a free outflow boundary 
condition is employed at the outer one. 



2.2. Input physics 

Self-gravity is taken into account by solving Poisson's equation 
in its integral form, using an expansion into spherical harmon- 
ics as in ,Muller & Steinmetz| ( p"995j . The monopole term of the 
po tential is corrected fo r gen eral relativistic eff'ects as described 
in IScheck et aT] ([2006^ and 'A rcones et aL] ( [2007] ). The inner 
core of the PNS with densities well above those of the neutri- 
nospheric layer is excised and replaced by a point mass at the 
coordinate origin. The cooling of the PNS is then described by 
neutrino emission properties (luminosities and mean spectral en- 
ergies) that are prescribed as boundary conditions at 7?ib. The 
contraction of the PNS is mimicked by a radial grid movement 
with defined velocity according to |Scheck"et a l. (2006); detailed 
information on the chosen parameters can be found in Sect. 2.1 



of the recent paper by [MuUer et al . | ( [20 1 2b| ) . There is no mass 
flow through the inner grid boundary as long as we follow the 
evolution of the PNS surface layers. 

"Ray-by-ray" neutrino transport and neutrino-matter interac- 
tions are approximated as in Scheck et al., (2006 ) by radial inte- 
gration of the one-dimensional (spherical), grey transport equa- 
tion for Vg, Vg, and heavy-lepton neutrinos and all angular grid 
directions (0, (p) independently. This approach allows us to take 
into account angular variations of the neutrino fluxes produced 
by the matter accreted onto the newly forming NS. The accre- 
tion luminosity adds to the neutrino fluxes imposed spherically 
symmetrically at the inner grid boundary. The neutrino spectra 
are assumed to have Fermi-Dirac shape with chemical potentials 
that are equal to the equilibrium values in optically thick regions 
and constant, limiting values in the free streaming regime. The 
spectral temperatures are computed from inverting the ratio of 
neutrino energy and number fluxes, whose transport is solved si- 
multaneously. Details can be found i n jS check et al. ( 2006|). T he 



tabulated equation of state (EoS) of |Janka & Muli er 
used to describe the stellar fluid 



|T9%) i 

It includes arbitrarily degen- 
erate and arbitrarily relativistic electrons and positrons, photons, 
and four predefined nuclear species (n, p, a, and a representative 
Fe-group nucleus) in nuclear statistical equilibrium. 



2.3. Nucleosynthesis 



In order to follow approximately the explosive nucleosynthe- 
sis, we solve a small a-chain reaction network, similar to the 
network described in "Kifoni dis et al.| ( |2003| ). For our simula- 
tions we choose a subset of nine a nucleiT^e, ^^C, ^^O, ^^Ne, 
^^Mg, 2^Si, 40Ca, "^^Ti, ^^Ni, and an additional "tracer nucleus". 
We omit ^^S, ^^Ar, ^^Cr, and ^^Fe from the complete a-series 
in order to reduce the network size and thus to save comput- 
ing time. The two gaps between ^^Si-'^^Ca and "^"^Ti-^^Ni are 
bridged by combining the intermediate reaction steps and es- 
timating the new eff'ective reaction rates by taking the slowest 
rate of the intermediate reactions. The tracer nucleus is pro- 
duced via the reaction ^^Ti(3Qf,7)^^Ni within grid cells whose 
electron fraction Ye is below 0.49. Such conditions are found in 
the neutrino-heated ejecta of our models. The tracer nucleus rep- 
resents iron-group species that are formed under conditions of 
neutron excess. It thus allows us to keep track of nucleosynthesis 
in neutron-rich regions. However, it should be noted that some 
fraction of the tracer material may actually be ^^Ni, because 
our approximations in the neutrino transport tend to underesti- 
mate Ye in the neutrino-heated ejecta, while more sophisticated 
energy-dependent neutrino transport yields (slightly) proton-rich 
conditions in the neutrino-processed material expelled during the 
early explosion ( e.g.,|Pru et et al. 2006; |Frohlich et al.|2006HFis-| 
Icher et al.|[20TQl |Muller et al.,2012a} . The network is solved 
in grid cells whose temperature is within the range of 10^ K - 
8 X 10^ K. We assume that all nuclei are photo-disintegrated to 
Qf-particles at temperatures above 8 x 10^ K. This is the high- 
temperature nuclear statistical equilibrium composition compat- 
ible with our a network in the absence of free neutrons and pro- 
tons. We neglect feedback from the detailed network composi- 
tion to the EoS and the hydrodynamic flow. This is an accept- 
able approximation because at SN conditions the contributions 
of nuclei to pressure, energy density, and entropy are dwarfed by 
those of electrons, positrons, and photons. 
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2.4. Long-time simulations 

To follow the propagation of the SN shock wave until very late 
times, we map our results after 1.1-1.4 s onto a new computa- 
tional grid, whose inner and outer radial boundaries are placed 
at 500 km and 3.3x10^ km, respectively. The latter value is near 
the stellar surface. At the inner grid boundary we assume the in- 
flow of a spherically symmetric neutrino-driven wind, which is 
considered to be the consequence of ongoing neutrino emission 
and corresponding energy deposition near the PNS. We do not 
directly simulate neutrino-matter interactions during our long- 
time runs, because details of the neutrino physics become less 
important during the late evolution phases when the explosion 
conditions are already determined. By including a wind inflow, 
however, we continue to account for the consequences of neu- 
trino heating at the PNS surface. Since the imposed wind has 
spherical symmetry it will not afl'ect the anisotropy of the lead- 
ing ejecta that we intend to investigate. The hydrodynamic in- 
flow quantities are obtained as the angular averages of the flow 
properties at 500 km at the time of the mapping and kept constant 
for another 2 s of evolution. 

The wind is replaced by a free outflow boundary condition 
afterwards. At the outer boundary we apply a free outflow con- 
dition at all times. Since we keep the angular grid resolution 
unchanged but increase the number of radial zones, our com- 
putational grid grows to 1200(r) x 41(6) x 137(0) x 2 zones. 
The radial resolution is approximately Ar/r ^ 0.01. We succes- 
sively discard the innermost radial zone and thus move the inner 
grid boundary outward to a new location whenever its radius be- 
comes less than 2% of the minimum shock radius. Moving the 
inner grid boundary to a larger radius helps relaxing the CFL 
timestep and allows us to follow the evolution to very late times 
with acceptable computational costs. 

Also during the long-time runs we consider nucleosynthe- 
sis processes and account for the self-gravity of the stellar gas 
by a Newtonian description. We switch to the tabulated EoS 
of |Timmes & Swesly] pOOO| ), appropriate for the lower density 
and temperature values pres ent in the outer stellar layers (not 
covered by the EoS table of |Janka & Muller||1996| ). Moreover, 
besides considering the contributions from arbitrarily degenerate 
and relativistic electrons and positrons and from a photon gas, it 
treats the nucleonic gas components by a mix of ideal Boltzmann 



gases with 11 nuclear species: p, "^He, ^^C, ^^O, ^^Ne, "^"^M^ 



^ai .p^ci^.. p, I2n 16o 20xT. 24, 

^^Si, ^^Ca, ^^Ti, ^^Ni, and X. This is suflftcient for a reasonably 
good representation of the composition in the outer core, mantle, 
and envelope of the progenitor. 



2.5. Computed models 

We investigate three 15M0 models and one 20 M© progenitor 
denoted as W15, L15, B15, and N20, respective ly. W15 is 
based on the nonrotating 15Mo model sl5s7b2 of |Woosley & 
Weaver (1995), LI 5 corresponds to a star evolved by [Limongi 



et al.,(,2000,), B15 is a blue supergiant (SN 1987A) progenitor of 
Woosley et al.| ([1988|), and N20 was c omputed up to core col- 
lapse by'Shigeyama & Nomoto|( [1990| . The 15 M© progenitors 
W and L were followed through collapse to 15 ms after bounce 
with the Prometheus -Vertex code in one dimension (A. Marek 
and R. Buras, private communication). Using the same code, the 
progenitor N20 was evolved until 1 1 ms after bounce (A. Marek, 
private communication), wherea s the model B 15 was collapsed 
in a ID simulation by S. Bruenn ( |Bruenn|1993 ). To break spher- 
ical symmetry, random seed perturbations with an amplitude of 
0.1% (and cell-to-cell variations) are imposed on the radial ve- 



locity (Vr) fielcj^ Explosions with chosen energy are initiated by 
neutrino heating at a rate that depends on suitable values of the 
neutrino luminosities assumed at the lower boundary. 

We computed in total 20 models varying the neutrino lumi- 
nosities imposed at the inner grid boundary and thus the result- 
ing explosion energies (Table [T]). Some of our models (W15-1 
and W15-2; B15-1 and B15-2) difl'er only by the initial seed per- 
turbations. The evolution including our treatment of neutrino 
transport is followed until 1.3 s after bounce for the W15 and 
N20 models, while the LI 5 simulations are continued until 1.4 s 
postbounce and the B 15 runs until about 1.1s after bounce. The 
time at which we stop our calculations difl'ers from progenitor 
to progenitor depending on how fast the explosion sets in. Ex- 
plosions in the L15 models are the most delayed ones, while the 
B15 models show the fastest explosion. The simulations with 
detailed neutrino treatment are stopped when the maximum SN 
shock's radius is close to 7?ob. Models which are simulated with 
the reduced angular grid resolution of 5° are indicated by the suf- 
fix "Ir" appended to their model names. In addition, a subset of 
9 models is continued into the long-time evolution as described 
in Sect. |2.4| in order to study the further acceleration of the PNS 
and the explosive nucleosynthesis behind the outgoing shock. 

3. Neutron star kicks 



Wongwathanarat et"aL] p010b') presented results of NS kicks for 
a small set of 3D computations conducted until 1.3-1.4 s after 
bounce. Here, we describe the evolution of these models for a 
much longer time as well as more results of a significantly en- 
larged model set that includes runs for three 15 M© progenitors 
and the 20 M© star. Moreover, we analyse our numerical results 
in the light of theoretical considerations and simple analytical es- 
timates (extending the brief discussion in Wongwathanarat et al. 
2010b|and sharpening ar guments made by Scheck et al.||2006" 
and |Nordhaus et al.|2010| ). We also evaluate our simulations for 
recoil velocities caused by anisotropic neutrino emission in com- 
parison to the far dominant kicks by asymmetric mass ejection. 

3.1. Simulation results 

Figure [2] presents snapshots of model W15-2 at four difl'erent 
epochs: 0.140, 0.248, 0.515, and 1.30 s after bounce. Violent 
buoyancy and mass overturn is triggered by neutrino heating 
around the NS. At roughly 100 ms, Rayleigh-Taylor mushrooms 
begin to grow from the imposed seed perturbations in the con- 
vectively unstable layer that builds up between the gain radius 
(that separates neutrino-cooling and -heating regions) and the 
stalled SN shock. These high-entropy bubbles evolve in a highly 
non-stationary manner, start rising, growing, merging, partially 
collapse again and emerge once more to inflate to larger sizes. 
Along with the slowly expanding shock their angular diameter 
increases, adjusting roughly to the radial extension of the con- 
vective layer. Aided by convective overturn and global shock 
motions, possibly associated with activity due to the SASI, the 
delayed, neutrino-driven explosion sets in at about 250 ms post- 
bounce. 

To identify the explosion epoch, we perform a time- 
dependent evaluation of the explosion energy £^exp defined as the 
sum of the total (i.e., internal plus kinetic plus gravitational) en- 
ergy of all grid cells where this energy is positive. We define the 



^ Unfortunately, 3D progenitor models with self-consistent informa- 
tion about the asymmetries present in the stellar core before collapse 
are not available yet. 
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140 ms 





Fig. 2. Entropy-isosurfaces (left) of the SN shock (grey) and of the high-entropy bubbles (green), and entropy distribution in a cross-sectional 
plane (right) at ^ = 140,248,515 ms, and 1.3 s for model W15-2. The viewing direction is normal to the xz-plane. Length scales are given by 
the yardsticks below each image. Small Rayleigh-Taylor mushrooms start growing at ~ 100 ms (top panels). The rising high-entropy bubbles 
merge and rearrange to form larger bubbles by the time the explosion sets in (second from top). The NS starts accelerating due to the asymmetric 
distribution of the ejecta, the acceleration reaching its maximum at ~500ms (third from top). At this epoch, the ejecta show a clear dipolar 
distribution with more dense, low-entropy material concentrated in the hemisphere containing the kick direction. The final NS kick direction has 
already found its orientation, and its projection onto the xz-plane is shown by the white arrows pointing to the lower right direction. The bottom 
panels show the entropy structure when the essentially spherically symmetric neutrino-driven wind has developed, visible as the green central 
region enclosed by the highly aspherical wind- termination shock (lower right panel). 



Article number, page 5 of 



25 



Table 1. Explosion and NS properties for all models (see Sect. 3.1 for definitions) 



Model 
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explosion time ^exp as the moment when the computed explosion 
energy exceeds 10^^ erg. This roughly coincides with the time 
when the average shock radius exceeds 400-500 km. 

This is also the time when the NS develops a growing kick 
velocity along a well defined direction. While the kick direc- 
tion was randomly fluctuating before, the acceleration vector 
now settles and the NS velocity grows (see Wongwathanarat| 
[eTal . 2010b). A cross-sectional view of the entropy distribu- 
tion clearly shows a dipolar asymmetry with dense, low-entropy 
material concentrated at small radii primarily in the hemisphere 
to which the kick vector is pointing (Fig. |2] third row and also 
Fig. [3] upper row). At a later stage, the pattern of the nonradial 
structures freezes and the ejecta continue to expand in a nearly 
self- similar manner. An entropy slice shows the emergence of an 
essentially spherically symmetric neutrino-driven wind enclosed 
by an aspherical wind-termination shock. On the side where the 
explosion is stronger, which manifests itself by faster shock ex- 
pansion and higher entropy behind the SN shock (see the 10-11 
o'clock position in the bottom right panel of Fig. [2] and the up- 
per right panel of Fig.|3]), the wind is shocked to higher entropies 
because it passes the termination shock at larger radii and thus 
with higher velocity (see also ,Arcones et al.|2007| and ,Arcones 
|& Janka|20lT] ). 

Table[T] summarizes explosion and NS properties for all com- 



puted models: Mn- 



'^exp? 



^exp? 



Qfky, Jns, o^sk, and 



Tspin denote the NS mass, time of the onset of the SN explosion, 
explosion energy, NS kick velocity by the gravitational tug-boat 
mechanism, corresponding NS acceleration, additional neutrino- 
induced kick velocity, angle between gravitational and neutrino- 
induced NS kick, NS angular momentum, angle between NS 
spin and kick directions, and estimated final NS spin period. 
While results for all models are provided after 1.1-1.4 s of post- 



bounce evolution, v^^^^ and a^^^ give the NS velocity and accel- 
eration at 3.1-3.4 s after bounce for the subset of the simulations 
that was continued until this later time. 

The NS mass Mns is defined as the baryonic mass enclosed 
by the NS radius 7?ns, which is considered as the angular depen- 
dent radial location where the density is 10^^ g cm"^. In estimat- 
ing the NS recoil velocity, we assume conservation of the linear 
momentum of the whole progenitor star. Since the high-density 
core of the PNS is excised in our simulations, the NS cannot 
move when it exchanges momentum with the surrounding mat- 
ter. Thus it behaves like an object with infinite inertial mass, 
absorbing momentum in analogy to a wall that reflects a bounc- 
ing ball. Nevertheless, following |Scheck e Tal. (2006), we can 
estimate the NS velocity Vns from the negative of the total linear 
momentum of the ejecta gas as 



Vns(0 = -Pgas(0/Mns(0, 



(1) 



where Pgas = C °^ pv is the integral of the gas momentum 
over the computational volume between PNS radius and outer 
grid boundar}]^ 

In our simulations we find NS kick velocities ranging from 
about 30kms-^ (L15-3) to roughly 440 km s"^ (W15-6) at 1.1- 
1 .4 s after bounce. The NS acceleration a^^ is obtained by numer- 
ically diflTerentiating the NS velocity calculated from Eq. ([T]) and 
listed in Table [T] Because of the nonvanishing acceleration, the 
recoil velocity at 1.1-1.4 s p.b. is usually not the final NS veloc- 
ity. Only models ending with a low kick (like LI 5-2 and LI 5-3) 
may have reached nearly the asymptotic value already after this 
period of time. In the three long-time runs W15-1, W15-2, and 



Note that there is no momentum flux across the outer grid boundary. 
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Fig. 3. Entropy-isosurfaces {left) of the SN shock (grey) and of the high-entropy bubbles (green), ray-casting images of the density {middle), and 
entropy distribution in a cross-sectional plane {right) for models W15-1 {top) and LI 5-2 {bottom) at about 1.3 s and 1.4 s after core bounce. The 
viewing directions are normal to the plane spanned by the NS kick and spin vectors of each model, which also defines the plane for the entropy 
slices. The SN shock has an average radius of 13,000 km and 14,000 km for models W15-1 and L15-2, respectively (a length of 5000 km is given 
by yardsticks below the left images). The kick and spin directions are indicated by the white and black arrows, respectively, in the middle figures. 
The NS locations are marked by black crosses in the right plots. The images in the middle correspond roughly to the projections of the density 
distribution on the observational planes. Dilute bubble regions are light-colored in white and yellow, while dense clumps appear more intense 
in reddish and bluish hues. The purple circular areas around the NS represent the dense inner region of the essentially spherically symmetric 
neutrino-driven wind. The wind is visible in green in the right images and is bounded by the aspherical wind-termination shock. The NS is clearly 
displaced from the geometrical center of the expanding shock in the direction of the kick vector pointing to the lower left for model W15-1. The 
NS in model L15-2 has a much smaller kick velocity than in model W15-1, and does not show any clear displacement but remains located roughly 
at the center of the expanding ejecta shell. 



W15-6, we find a significant increase of the recoil velocity be- 
tween 1.3 and 3.3 s, and even then the acceleration continues, 
though at a considerably smaller rate (compare the values of Vns 
and ^ns with and without the superscript "long" in Table [T]). Al- 
though we evolved several of the long-time models to much later 
times we do not compute the NS velocities beyond 3.4 s, because 
the subsequent shift of the inner grid boundary and the applica- 
tion of a free outflow condition, can lead to anisotropic momen- 
tum flow across the boundary can occur. The latter would have to 
be accounted for with high precision if the argument of momen- 
tum conservation were to be used for estimating the NS recoil 
also later on. 



In Fig. [3] the asymmetry of the SN shock, the structure of 
the high-entropy bubbles, and the inhomogeneous, dense ejecta 
shell behind the expanding shock are compared for two models 
at 1.3-1.4 s. One model, W15-1 (top), has a large kick, the other 
model, LI 5-2 (bottom), a small one. The analogous combination 
of plots was provided for model W15-2, our second highest kick 
case, in Fig. 3 of , Wongwathanarat et al.| ( |2010bj . The three pan- 
els display diflTerent aspects of the anisotropic mass distribution. 
While in the left image the geometry is visible in overview by 
three-dimensional entropy-isosurfaces, the middle panel gives 
an impression of the asymmetry of the gas (column) density vi- 
sualized by a ray-casting technique. The viewing direction is 
perpendicular to the plane defined by the NS kick (white) and 
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time [s] 

Fig. 4. Long-time evolution of the NS velocity (top) and of the corre- 
sponding NS displacement (bottom) for models W15-6, W15-1, W15-2, 
L15-1, and L15-2. The displacement of the kicked NS from the coor- 
dinate center, although not taken into account in our simulations, can 
be estimated by time integration of the NS velocity. Within 3.3 seconds 
the NS in model W15-6 exceeds a velocity of 700 km s~^ and continues 
to be accelerated at a modest rate. Note also that the NS velocities in 
models LI 5-1 and LI 5-2 have nearly reached their terminal values after 
3 seconds. In all cases the NS shift relative to the coordinate center is 
negligibly small before the explosion takes off, and even in model W15- 
6 with its most extreme NS kick it is always dwarfed by the diameter 
of the shock wave and postshock ejecta, which expand with typically 
~10,000kms-^ 



spin (black) vectors. The right plot shows the entropy distribu- 
tion in that plane. The ray-casting images visualize the density 
projected along the lines of sight. They show that the kick vec- 
tor points to the region with the highest concentration of dense 
clumps. Moreover, the NS in model W15-1 is visibly displaced 
from the geometrical center of the ejecta and of the SN shock 
in the recoil direction; this effect is even more pronounced for 
W15-2. In contrast, in the low-kick case of L15-2 the high- 
entropy bubbles are of relatively similar sizes on all sides. The 
expansion of the SN shock is therefore nearly equally strong in 
all directions, i.e. the shock has no strong dipolar deformation 
(though a quadrupolar asymmetry mode is visible), and the NS is 
close to the center of the ejecta shell. Correspondingly, the wind- 
termination shock is relatively spherical in this model. Note that 
there are denser, clumpy structures near the one o'clock posi- 
tion of model LI 5-2 but also similarly massive, though less con- 
centrated filaments in the opposite hemisphere. The explosion 
asymmetry has a dipolar component that is not large enough to 
produce a high pulsar kick. 

Figure |4] displays the growth of the NS velocities in five 
of our nine long-time runs over the whole period of 3.3-3.4 s. 
In the most extreme cases, models W15-6 and W15-2, the NS 
velocity reaches an asymptotic value over 700 km s"^ and near 



600 km s" , respectively. Several of the W15-models therefore 
develop kicks clearly higher than the peak of the measured ve- 
locity distribution of young pulsars. Only more simulations for 
a wider variety of progenitors, however, will be able to provide 
answers to the questions whether velocities of even 1000 km s"^ 
and more can be obtained by the discussed gravitational tug-boat 
mechanism, and whether and how such extreme recoils might be 
connected to particular features of the progenitor stars or espe- 
cially favorable conditions of the explosion itself. Analytic es- 
timates in Sect. |3.4| will show that the mechanism is likely to 
be viable of producing kicks well in excess of 1000 km s"^ for 
reasonable assumptions about the explosion asymmetry. Yet it 
is not clear that the magnitude of the kick is merely a stochas- 
tic result of the chaotic growth of asymmetries by hydrodynamic 
instabilities in the SN core. Some trends of our results seem to 
indicate that certain properties of the progenitor and the explo- 
sion could be supportive for large kicks (see also Sect. |3.5| ). 

3.2. Mode analysis 

The causal connection between the explosion asymmetries, es- 
pecially the density asymmetries in the ejecta shell, and the natal 
NS kicks, which is visually suggested by Figs. [2] and |3] as well 
as by Fig. 3 of |Wongwathanarat et al.| ( [2010b| ), can be quantita- 
tively supported by an expansion of variables that characterize 
the mass distribution in the postshock layer into spherical har- 
monics. 

Testing different possibilities, we ascertained that the most 
suitable quantity for this purpose is the "surface mass density" 

D(6,(P)^ drp(r), (2) 

where Rs(0, (p) is the direction-dependent radius of the SN shock. 
The integrand of Eq. ^ means that mass elements Am at differ- 
ent distances from the NS are weighted according to Am/r^, re- 
flecting the importance of their contribution to the NS accelera- 
tion by gravitational effects. In contrast, the shock radius Rs(0, (p) 
or the integrated mass per solid angle, 

W,(f>)= I drrV(r), (3) 

put too much weight on the region just behind the shock, which 
is not necessarily the most relevant one for accelerating the 
NS. The NS kick is mainly determined by the density inhomo- 
geneities in the deeper ejecta regions. Although a large kick 
velocity is typically also linked to a pronounced £ = I displace- 
ment of the shock surface relative to the grid center (in our most 
extreme cases the dipole amplitude can reach nearly 10% of the 
monopole one, see Fig. |5j left panel), we find that the inverse 
conclusion does not always hold: A dominant £ = 1 shock as- 
phericity does not necessarily imply a strong NS recoil. In fact, 
we observe the £ = 1 spherical harmonics of the shock surface 
to be larger than all higher-order multipole amplitudes at late 
postbounce times for the majority of our models (except for the 
few cases like L15-2, L15-3 or N20-2, all of which are among 
our runs with the lowest NS kicks). For a discriminating diag- 
nostic quantity we therefore prefer to evaluate D(6, 0), whose 
mode amplitudes are more conclusive about ejecta asymmetries 
that are causal for large recoil velocities of the compact remnant. 
Large kicks correlate with an ^ = 1 mode of D(0, (p) that clearly 
dominates the higher modes during (at least) a longer period of 
the postbounce evolution, whereas for cases with small kicks the 
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Fig. 5. Time evolution of the normalized spherical harmonic expansion coefficients for diff'erent {i, m) modes computed for the shock surface 
Rs{6, 0) {left) and for the surface mass density D{6, 0) (Eq.[2j right) of models W15-2 (top) and L15-2 (bottom). The polar axis of the coordinate 
system chosen for the expansion is aligned with the final NS kick direction. Note the dominance of the {i, m) - (1,0) mode for model W15-2 with 
the highest kick velocity, while this dominance is absent in the case of model LI 5-2. 



1=1 mode is not much stronger at any time than the other 
modes. 

The decomposition of D{6, (p) in spherical harmonics is writ- 
ten as 



CO i 



(4) 



^=0 m^-e 



where are the expansion coefficients and the spherical har- 
monics are functions of the associated Legendre polynomials 
as 

Ff(6>,0) = ^fP^(cos6/)^^'^'^, 
with 



(5) 



(6) 



Multiplying Eq. (|4]) by the complex conjugate of the spherical 
harmonic, 7^*, and integrating over the solid angle, the expan- 
sion coefficient is found to be 



Jo Jo 



de smeD(0,(f))Y^\e,(f)). 



(7) 



We choose the coordinate system for the expansion such that 
the polar axis is aligned with the final NS kick direction. Be- 
fore performing the integration of Eq. (|7]), we interpolate, us- 
ing bi-linear interpolation, the surface mass density computed 



on the Yin- Yang grid onto the new coordinate system. Figure [5] 
right panel displays the resulting expansion coefficients of dif- 
ferent (£, m) modes versus time for models W15-2 and L15-2. In 
the former, high-kick model the amplitude of the (^, m) = (1, 0) 
mode reaches its maximum at about 400 ms after bounce, i.e., 
approximately at the same time when the NS acceleration be- 
comes maximal (Fig.|6]). After this time, at ^ ^ 0.4 s postbounce, 
the ii,m) = (1, 0) mode clearly dominates. Such a behaviour is 
absent in the low-kick model LI 5-2. 



Consistent with what is visible in Figs. |2] [3] and Fig. 3 of 
[Wongwath anarat et al.| ([2010b), the mode analysis shows that a 
dominant f = 1 (dipolar) asymmetry of the mass distribution in 
the ejecta correlates with a high pulsar kick along the dipole axis. 
In model LI 5-2 all of the £ = 1 and £ = 2 modes displayed in the 
lower right panel of Fig. [5] exhibit similar amplitudes and thus 
the gravitational pull of the ejecta does not develop a preferred 
direction. The acceleration of the NS remains correspondingly 
low. This result provides strong support to the described grav- 
itational tug-boat mechanism for NS kicks, which is linked to 
the anisotropic ejection of matter in the first second of the super- 
nova explosion and which is mediated to the NS on significantly 
longer timescales by the asymmetric gravitational pull of the ex- 
panding ejecta gas. In the next section we will describe in more 
detail how the asymmetries are established and how the acceler- 
ation mechanism works. Moreover, we will present simple toy 
models to approximately estimate its possible strength in depen- 
dence of the properties of the supernova ejecta. 
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Fig. 6. Left: Time evolution of the NS acceleration in model W15-2. The acceleration computed by the assumption of momentum conservation, 
i.e., by the negative derivative of the total linear momentum of the gas surrounding the NS (red; denoted by "mom.consv."), is compared with 
the value obtained by the integration of all hydrodynamical and gravitational forces that act on a sphere of radius 1 .3/?ns around the PNS (black; 
denoted as "total"). Also shown are the absolute values of the individual contributions of all terms in Eq. (jSj, namely the momentum transfers by 
downflows (light blue), outflows (green), anisotropic pressure forces (purple), and the gravitational attraction by the asymmetric ejecta (blue). The 
last contribution clearly dominate the evolution over long periods. Right: NS velocity components in all coordinate directions and total value of 
the NS velocity. Note that the two lines in each panel, which correspond to our two methods of computing the NS recoil acceleration in the left 
panel, can hardly be distinguished from each other. 



3.3. Kick mechanism: theory and tests 

The NS acceleration that is associated with the asymmetric 
ejecta distribution developing during the onset of the SN explo- 
sion is achieved mainly by gravitational forces in combination 
with hydrodynamical forces. This will be detailed below. The 
acceleration proceeds in three steps: 

(1) When violent convective mass flows and possibly SASI 
sloshing motions conspire to stir the postshock layer, an 
anisotropy of the mass-energy distribution around the PNS is 
created. Convective downdrafts, channelling accretion flows 
from the stalled shock through the neutrino-heating region 
to the vicinity of the gain radius, get deflected to feed an 
asymmetric pattern of high-entropy bubbles. The bubbles 
form, grow, contract again, partly collapse, and reappear in 
a quasi-chaotic way to become smaller or larger, absorbing 
less or more neutrino energy. This stochastic bubble forma- 
tion, however, does not cause an appreciable recoil of the NS 
(see Fig . [6| and | Wong wathanarat et al . |2Q 1 Ob ) . 

(2) When the explosion sets in, the shock and postshock gas be- 
gin to expand aspherically, driven by the asymmetric infla- 
tion of the bubbles. The ejecta gas therefore gains radial 
momentum and its c.o.m. moves away from the coordinate 
origin (Figs. |2] [3] right panels): The ejecta shell acquires a 
net linear momentum because of difl'erent strengths of the ex- 
plosion in difl'erent directions. The initial energy and mass 
asymmetry is thus transformed into a momentum asymme- 
try by the conversion of thermal to kinetic energy through 
hydrodynamical forces. When the expansion timescale be- 
comes shorter than the timescale of lateral mixing, the asym- 
metric ejecta structures essentially freeze out. 

(3) Because of linear momentum conservation, the NS must re- 
ceive the negative of the total momentum of the anisotropi- 
cally expanding ejecta mass. Hydrodynamic pressure forces 
alone cannot achieve the NS acceleration (Sch eck et al.| 
[2006). As long as accretion downdrafts reach the NS, mo- 
mentum is transferred by asymmetric gas flows. Stronger 
accretion on the weaker side of the blast and more mass 



loss in the neutrino-driven wind on the other side cause a re- 
coil opposite to the main explosion direction. However, the 
largest kick contribution, which continues even after accre- 
tion has ended and after the wind outflow has become spher- 
ical (Fig. [2] bottom right panel, and Fig. [3] right panels), 
results from the gravitational pull of the anisotropic shell of 
ejec ta ([Scheck et al. 2006; Wongwathanarat et al.|2010b| and 
also |Nordhaus et al.|2010| ). 



Why can our models produce kick velocities in excess of 
~500kms"^ by hydrodynamically created explosion asymme- 
tries although such a mechanism appeared to be highly disfa- 
vored even for the extreme case of dipolar ejecta deformation in 
an analysis of 2D explosion models by| Janka & Muller| ( |1994| )? 
It should be noted that the kick scenario discussed in the lat- 
ter paper is fundamentally difl'erent from the one described here. 
|Janka & Muller| ( |1994| ) assumed the kicks to originate from an 
impulsive event, happening during a short period of a few hun- 
dred milliseconds around the onset of the anisotropic explosion. 
Correspondingly, in their estimates they considered the momen- 
tum asymmetry of the ejecta shell during this early phase, ex- 
pecting the kick to be established before the accretion is over. 
Indeed, also in the simulations described in the present paper 
the NS velocity at the end of the accretion phase is much lower 
than 500 km s"^ in all cases. However, in the scenario discussed 
here the NS acceleration is a long-duration phenomenon, which 
continues for several seconds even after accretion has ended, be- 
cause it is promoted by the long-range gravitational forces of the 
asymmetric, dense ejecta shell that expands away from the NS 
behind the outgoing SN shock. 



To elaborate on this point and to check our esti mates of the 
NS recoil velocity on the basis of Eq. ([T]), we follow [SchecketaL 
(|2006 ) and perform an analysis of the difl'erent contributions to 
the NS kick by summing up all forces a cting on a central sphere 
of radius ro = 1.37?ns (Nordhau s et al.|[2Q1 0, adopted a similar 
analysis for their recent 2D results). According to Eq. (5) of 
[Scheck et al.| (2006) the time-derivative of the NS momentum 
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can be deduced from the Euler equations as 



Pns^-Cp rdS-(p pyVrdS^l GMns-dm. (8) 

The corresponding NS acceleration is ans = Pns/^ns and the 
NS velocity Vns(0 = ^ns(t')- The first term on the rhs of 

Eq. ([5]) accounts for the variation of the pressure force around the 
sphere. The second term is the flux of momentum going through 
the sphere. This term can be decomposed into downflow and 
outflow contributions depending on the sign of the radial velocity 
Vr being negative or positive, respectively. The third term is the 
contribution from the long-range gravitational forces exerted on 
the NS by the anisotropic mass distribution outside the sphere. 
Note that a change of the radius of evaluation, tq, in a reasonable 
range around 1.3/?ns (i-^., in the neighborhood of the NS surface) 
leads to di fl^erent weights of the th ree summed terms (as noticed 
recently by |Nordhaus et al. |20 1 0} but still yields the same result 
for the integrated total eff'ect. 

Figure [5] compares a^s and Vns for model W15-2 as calcu- 
lated by the two difl'erent methods, i.e. from the requirement of 
total momentum conservation (Eq. [T]) and from the integration 
of the responsible hydrodynamical and gravitational forces ac- 
cording to Eq. ([8]). NS acceleration as well as velocity exhibit 
excellent agreement of the results obtained in the two difl'er- 
ent ways. Between -130 and -600 ms after bounce the forces 
by the anisotropic pressure distribution and by the gravitational 
pull both contribute to the total acceleration on significant lev- 
els, while downflows and outflows transfer sizable amounts of 
momentum mostly during transient phases of high activity when 
strong accretion streams impact the nascent NS anisotropically. 
At later stages, i.e., at ^ ^ 400ms when the asymmetry of the 
ejecta has been established, and in particular at ^ ^ 600 ms when 
accretion to the PNS has ceased and the neutrino-driven wind 
has become essentially spherically symmetric (Fig. [2]), the grav- 
itational force term starts to dominate and therefore accounts for 
the long-term increase of the NS kick. A major fraction of the 
final NS velocity is thus caused by the gravitational attraction of 
the anisotropic ejecta. 

Our analysis, compared to the 2D results of [Scheck et al.| 
2006| ), reveals that in the 3D case the direct hydrodynamical 



momentum transfer by gas downflows to and outflows from the 
NS is much less important than in 2D. While in 2D the toroidal 
nature of all structures favors the formation of large hemispheric 
difl'erences with inflow to the NS dominating in one hemisphere 
and outflow in the other (corresponding to an ^ = 1 geome- 
try), the 3D downflow and outflow pattern is more often and for 
longer periods of time characterized by higher multipoles. The 
momentum transferred to the NS by inflows is thus almost bal- 
anced by the momentum carried away from the NS by outflows 
on smaller angular scales, resulting only in a small contribution 
to the NS kick velocity. The gravitational tug by the asymmet- 
ric ejecta distribution is therefore clearly the dominant mediator 
of the NS acceleration in a time-integrated sense, although it is 
not necessarily the only important efl'ect, in particular not dur- 
ing the very early phase of the NS recoil and not in 2D. Sizable 
directional difl'erences of the explosion strength can occur and 
can lead to large ejecta inhomogeneities and asymmetries also 
in 3D (most pronounced in models W15-6 as well as W 15-1 and 
W15-2, see Figs. |2] and [3]). Large recoil velocities can therefore 
be produced despite the fact that the pattern of convective bub- 
bles and downdrafts seems to be characterized by the presence 
of higher-order multipoles. 

Note that allowing the PNS to move (instead of being at- 
tached to the coordinate center because of the excision of the 



central core in our simulations) is unlikely to change fundamen- 
tal aspects of our results. This expectation is based on a simple 
fact visualized in Figs.|4](bottom panel) andjTJ Even for the high- 
est computed kick velocities the displacement of the NS from 
the coordinate origin during the computed evolution of several 
seconds (at most -1500 km) is very small compared to the radial 
extension of the explosion shock and asymmetric ejecta distribu- 
tion (typically several 10,000 km). It is therefore not astonishing 
that our expectation was recently confirmed by [Nordhaus et al.| 
(|2010|[2012j) with 2D simulations, in which the NS was able to 
move. The minor influence of the NS motion had been tested 
before by Sche ck et al.| (^2006 ), who performed a large set of 2D 
simulations in which some of the kicks were similarly strong as 
the 3D cases discussed here. For doing the tests, [Scheck et al. 
(2006]) applied a Galilei transformation with velocity -Vns to the 
gas outside of the PNS (see Sects. 2.3, 7.2, Appendix B, and 
Table A.5 in the cited paper), thus mimicing the efl'ects of an 
accelerating NS moving through the surrounding medium. Nat- 
urally, because of the stochastic and chaotic growth of local and 
global asymmetries, the simulations developed values of Vns that 
could be difl'erent for a fixed progenitor and the same choices of 
boundary conditions (core neutrino luminosities and PNS con- 
traction) and initial perturbations. Nevertheless, the whole set of 
models with the NS motion taken into account did not exhibit 
any systematic diff'erences compared to the standard set with- 
out the Galilei transformation applied (cf. Tables A.l and A.5 in 
[Scheck et al.|2006| ). 



3.4. Analytic estimates of tlie kick magnitude 

What are the maximum kick velocities that can be mediated to 
the NS by gravitational forces of asymmetrically ejected matter? 
The maximum velocity of more than 700 km s"^ obtained in our 
present set of 3D models is no upper limit, and recoil veloci- 
ties beyond 1000 km s"^ were already found in 2D simulations 
by [Scheck et al.| ( [2Q06| ). Simple considerations can yield order- 
of-magnitude estimates that confirm this possibility. In the fol- 
lowing we discuss two diflferent, idealized cases that allow one 
to integrate the gravitational pull over the whole period of accel- 
eration (see also IWong wathanarat et al.|2010b| and alternative, 
simple estimates by Nor dhaus et al.|2010^ 



In the first case we consider a hemispheric asymmetry of the 
mass distribution in an expanding, spherical ejecta shell. The 
shell with an average radius is assumed to contain an extra 
clump-like mass Am in one hemisphere, while an equal mass 
deficit exists in the other hemisphere. A mass asymmetry of this 
kind may be established by asymmetric accretion flows during 
the violently convective postbounce phase. Downflows may, for 
example, be stopped and blown out again before they can reach 
down to the PNS surface, or only some part of the downflows 
is accreted onto the PNS while the rest is deflected and ejected 
again after neutrino heating. Let us assume that the shell with its 
time-dependent radius r^it) expands with a constant velocity Vs, 
starting at an initial radius If the NS is displaced by a distance 
s from the center of the shell (s = at ^ = 0), the gravitational 
forces lead to a NS velocity Vns = s via an acceleration a^s = s 
given by 



.. dvns ^ . 
s = — — = GAm 
dt 



1 



1 



(rs - s)^ (r, -h s)^ 



(9) 
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Fig. 7. Ejecta morphology of the high-kick models W15-2 (left) and W15-6 (right) at approximately 3.3 s after bounce. The NS (indicated 
by black spheres, enlarged by a factor of 20 for better visibility) are accelerated by the gravitational tug of the asymmetrically distributed, dense 
ejecta to velocities of 575 km s~^ and 704 km s~\ respectively (green arrows pointing to the right). The ejecta "clumps" with the highest densities, 
which are mainly responsible for the gravitational attraction, are visualized by the red-colored surfaces of constant density of lO^gcm"^. The 
small black dot left of the NS in both images marks the coordinate origin. The NS displacement relative to this location was computed by time 
integration of the NS velocity. After ~3.3 s it is only 1260 km in W15-2 and 1440 km in W15-6. The neutrino-heated, asymmetrically distributed 
inner ejecta are bounded by the semi-transparent, beige surface from outside and the semi-transparent, light blue surface from inside. The former 
surface is defined by a constant ^^Ni mass fraction of 5%, the latter surface conincides with the termination shock of the neutrino-driven wind. The 
outermost, semi-transparent green surface (cut open towards the observer for better visibility of the inner volume) corresponds to the SN shock 
wave. The average shock radius is approximately 36,000 km at the displayed time. A yardstick of 10,000 km along with a tripod indicating the 
coordinate directions are provided in the lower right corner of each plot. 



Using Ts = Ti + and assuming s ^ at all times, the integra- 
tion of Eq. ^ from t = Otot = oo yields: 



2GAm 



540 



km 

s 



Am_3 

n,7 Vs,5000 



(10) 



where Am is normalized by 10"^ M©, ri by 10^ cm, and Vg by 
5000 km s"^ A 10"^ M© hemispheric asymmetry of the shell ex- 
panding with a constant velocity of 5000 km s"^ can thus pull 
the NS to a velocity of 540kms~^ (this result was already men- 
tioned in Wongwathanarat et al. 2010b ). Ejecta asymmetries can 
therefore very effectively mediate a long-lasting acceleration of 
the NS, which is pulled in the direction of the higher mass con- 
centration. According to Eq. ([TO]) the kick becomes larger for a 
lower expansion velocity, because the gravitational tug from the 
asymmetric ejecta shell remains strong for a longer time. While 
a value of 5000 km s"^ tends to be on the high side for the post- 
shock ejecta shell during the relevant times, a clumpy anisotropy 
of 10"^ Mq means an asymmetry of the mass distribution of only 
in a shell that typically contains 10"^ Mq or more. Asym- 
metries of several percent appear well within reach and there- 
fore velocities of more than 1000 km s"^ as well. Note that the 
first-order terms of s/r^ cancel in Eq. ([9]). This means that the 
NS recoil is aflTected by the NS motion and displacement from 
the explosion center only on the level of second-order correction 
terms of s/r^ <$c 1. 

In the second case we consider a diflTerence of the expan- 
sion velocity of two clumps Am of the ejecta in both hemi- 
spheres instead of a hemispheric mass diflTerence. This means 
that we assume that these clumps have been accelerated diflTer- 
ently strongly and propagate away from the center of the blast 
according to rgj = n -\- v^^it and rs,2 = n -\- Vs,2^ with Vs,i < Vs,2. 



Taking s to be again the displacement of the NS from the blast 
center, the compact remnant in this situation experiences a grav- 
itational acceleration 



s = — — = GAm 
dt 



1 



1 



(^8,1 - s)^ (rs,2 + s)^ 



(11) 



To lowest order in s/r^ <^ 1, time integration from to oo leads 
to 



GAm Vs 2 - V| 



s,l 



n Vs,l Vs,2 



GAm Avs 
[km 

260 



Am_ 



n,7 Vs,5000 



Avc 

^ , (12) 



where we have introduced the definitions Avs = Vs,2 - Vs,i and 
4Vs,2. In the last expression we have again normalized 
Ti by 10^ cm, and Vs by 5000 km s"^ The accel- 



Vs = 

Am by 10"^ M©, 
eration is opposite to the faster expanding hemisphere, i.e. it is 
against the direction of the stronger explosion, and of the same 
order of magnitude as the result of Eq. ([TO]), if we assume a ve- 
locity asymmetry of 100%, Avs/Vs ~ 1- This, however, is on the 
very extreme side, and Avg/Vs < 0-5 seems more realistic. Ac- 
cording to Eqs. (fTTl) and ([12]) a velocity asymmetry is therefore 



less eflftcient in accelerating the NS than the mass asymmetry in 
Eqs. ([9]) and ([TO]). The first-order correction term in s/r^ <^ 1, 
which we have suppressed when going from Eq. (fTTl) to Eq. ([12]), 



reads dv^sV^^ = 2GAm s{r \ + r It is not necessarily small 



compared to the leading term dv^^Vd^ = GAms(r~f - r'^), but 
the corresponding acceleration amplifies the effect of the zeroth- 
order term, because the displacement by the kick brings the NS 
closer to the attracting slower parts of the ejecta. 
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Fig. 8. NS kick velocities and acceleration values (indicated by the 
length of the arrows attached to the symbols) versus explosion ener- 
gies for all computed models after 1.1-1.4 s of simulated postbounce 
evolution. Models corresponding to different progenitors are shown by 
different symbols as defined in the figure. 



A closer inspection of Eqs. ([TO]) and ([12]) clarifies the mean- 
ing and implications of the non-impulsive nature of the described 
gravitational tug-boat mechanism for NS acceleration. The in- 
stantaneous linear momentum carried by asymmetrically dis- 
tributed clumps with representative masses and velocities, Pgas ~ 
2AmVs ~ 2 X 10^^ Am_3Vs,5ooo gcms"^ is one to two orders of 
magnitude smaller than the momentum associated with the esti- 
mated possible NS kick, Pns = ^ns,i.5 Vns,8 ~ 3 X 10^^ gcms"^ 
with Mns,i.5 = Mns/1.5 M©. Efficient NS acceleration can there- 
fore not be explained by considering just the instantaneous mo- 
mentum of the clumps at any time. Exactly this is the problem 
connected with purely hydrodynamically mediated kicks where 
the NS receives its recoil by hydrodynamical forces during the 
relatively short (hundreds of milliseconds) phase of anisotropic 
mass accretion and ejection. Correspondingly, Janka & Muller 
1994| ) showed that even with the most extreme assumptions for 



the explosion asymmetry these more impulsive hydrodynami- 
cal kicks are unlikely to lead to recoil velocities of more than a 
few 100 km s"^ Instead, stronger kicks require a long-lasting 
accelerating force, which is provided by the gravitational in- 
teraction between the NS and the asymmetric ejecta. For this 
gravitational tug-boat mechanism to work, it is necessary to sus- 
tain the momentum of the dense clumps in the ejecta shell over 
many seconds, despite the gravitational transfer of clump mo- 
mentum to the NS. This requires a continuous clump acceler- 
ation which can happen by hydrodynamical forces within the 
inhomogeneous ejecta, whereby internal energy of the ejecta is 
converted into kinetic energy provided the ejecta do not expand 
ballistically yet. The latter requirement limits the radial scales 
and time intervals that can contribute to the pulsar kick. The 
clumps need to be hydrodynamically coupled to the energetic 
environment they are embedded in. This guarantees that they are 
continuously accelerated while they transfer momentum and en- 
ergy to the NS through their long-range gravitational attraction. 
The energy thus pumped into the motion of the NS can become 
quite significant, ^k,ns = l^nsVns ~ lO'^^Mnsj.sv^^ 8 erg. This, 
however, is still a small fraction of the ~10^^ erg of a normal SN 
explosion. 



3.5. Progenitor dependence of neutron-star kicks 

The data of Table [T] seem to suggest an interesting trend in 
dependence of the progenitor models: All computed cases for 
the W15-progenitor yield NS recoil velocities that clearly ex- 
ceed 200km/s, and in all cases the acceleration is still relatively 
high at the end of the simulation. Also two of the L15 mod- 
els (the low-resolution run L15-4-lr and LI 5-5) develop kicks of 
~200km/s and higher with still large acceleration at the end of 
the simulated evolution, but all models of the B and N progeni- 
tors produce clearly slower NSs. 

Both model subsets of strong and weak kicks exhibit a range 
of explosion energies without any clear correlations between 
energy and kick velocity (Fig. [S]). Note that for each investi- 
gated progenitor the number of models is still limited and in- 
cludes either a wider range of kick velocities for high as well 
as low explosion energies (W15, LI 5) or little variation of the 
kick velocities for all tested energies (B15, N20). We refrain 
from claiming a lack of high-kick cases for strong explosions 
(£exp ~ 1.2 X 10^^ erg) because no W15 simulation was done in 
this energy range and a trend could be connected with progeni- 
tor dependencies. In general, our present results are compatible 
with conclusions drawn by |Scheck et al.| (f2006 ) on the basis of 
a large set of more than 70 simulations in two dimensions and 
for three diflferent progenitor models with and without rotation. 
Scheck et"aL] ( [20061 ) observed kick velocities on the low and high 
sides of the velocity distribution for all investigated progenitors 
and for less as well as more energetic explosions. 

Presently we have no clue which model conditions, i.e., 
which progenitor or explosion or NS properties, could be fa- 
vorable or disfavorable for high kick velocities. Neither are 
W15-models with their typically higher kicks distinct concern- 
ing explosion energy or assumed NS contraction, nor is their 
progenitor structure clearly special in any aspect compared to 
the other stars. For example, L-models have the highest density 
and most shallow density profile outside of the iron core and thus 
the largest mass accretion rates after core bounce with the slow- 
est decay in time. Despite their longest delays until the onset of 
the explosion, however, they do not stick out concerning their 
NS kicks. 

The determination of clear trends in dependence of progen- 
itor, explosion properties, and remnant properties will therefore 
require a more extended investigation with a systematic variation 
of the involved degrees of freedom. This will need many more 
3D models than currently available, which is computationally 
extremely demanding, even with the approximations of the neu- 
trino transport used in our work. At the present time we therefore 
refrain from drawing premature conclusions or getting entrapped 
in unfounded speculations. 



3.6. Contribution by anisotropic neutrino emission 

A directional asymmetry of the neutrino flux radiated by the 
nascent NS leads to anisotropic energy and momentum loss from 
the compact remnant. This can provide a kick of the compact 
remnant, too. The contribution of this neutrino-induced kick to 
the NS velocity, Vns,y, is equal to the negative of the momentum 
carried away by the neutrinos divided by the NS mass. 



(13) 



This recoil can be expressed in terms of the total neutrino energy 
loss, Ey, by introducing dimensionless parameters kyj that de- 
scribe the asymmetry of the neutrino momentum loss in all three 
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Fig. 9. Neutrino flux asymmetry at ~150, 300, 600 ms after bounce 
for model W15-3 as representative case. The 4:7r-maps show the rel- 
ative variations AF/(F) of the total (i.e., sum for y and v of all three 
neutrino flavors) neutrino energy flux in all directions. Hot spots (yel- 
low and red) during the postbounce accretion phase (t < 400 ms; top 
and middle panels) are connected to strong convective downflows in 
the neutrino-heated postshock layer. The typical scale of the corre- 
sponding high-emission areas is 20-30 degrees, and the amplitude of 
variations ranges around 10-30% during most of the time, though rare, 
single maxima can reach up to ~70%. After the accretion has ended 
(t ~ 400-500 ms), the flux variation pattern is determined by the cell 
structure of the convection inside the proto-neutron star (i.e., below the 
neutrino sphere) and the variation amplitute decreases to a few percent 
(bottom panel). The number in the lower left corner of the plots (dS/dt) 
gives the total neutrino energy loss rate as flux integral over the whole 
emitting surface. 



coordinate directions: 



(14) 

where c is the speed of light. The vector components kyj are 
computed as 

dtj)dSnjPij'i-J^ dtj)dSh¥J . (15) 

Here Pij are the components of the neutrino pressure tensor P; 
Pij defines the flux component j of the neutrino momentum com- 



ponent /. Fy is the neutrino energy flux, d^" the surface element 
on the radiating (neutrino)sphere, n (with components nj) the 
unit vector perpendicular to the surface, and the time and sur- 
face integral in the denominator yields Ey. From Eq. ([14]) an 
order-of-magnitude estimate of the NS velocity is derived as 



3.3 X 10^ ky 



1.5 Mc. 



3 X 1053 erg M^, 



cm 

s 



(16) 



This yields a recoil velocity of Vns,y ~ 330 km s"^ for ky = 0.01. 

In principle, ultrastrong magnetic fields could lead to large- 
scale, potentially strong dipolar asymmetries of the neutrino 
transport out of the hot PNS, e.g. by changing the neutrino- 
matter interaction opacities, coupling to a neutrino magnetic mo- 
ment, or aff'ecting neutrino-flavor conversions via the MSW ef- 
fect (e.g., Bisnovatyi-Kogan 19961 [ Arras & Lai| 1999t|Kusenko| 
|2009l |Kishimoto||201 li niowever, little is known about the 
strength and structure of the magnetic field inside a nascent NS. 
It could well be chaotic due to convective activity, in which case 
a large global asymmetry with dominant low-order multipoles is 
unlikely. Not counting on support by strong dipolar B -field com- 
ponents, a directional bias of the neutrino emission of even only 
one percent is extremely difficult to obtain. 

Asymmetric neutrino emission can be caused by convective 
overturn in the postshock layer and by convective activity inside 
the NS, i.e., below the neutrinospheric region. In the postbounce 
accretion phase, asymmetric neutrino emission is dominated by 
low-entropy accretion in the neutrino-heating layer between gain 
radius and shock. Hot spots of high neutrino emission are con- 
nected to the strongest accretion funnels. They encompass typ- 
ical angular scales of 20-30 degrees and are scattered over the 
whole emitting surface (Fig. |9] top panel). The low-emission 
regions, which surround the hot spots, trace the distribution of 
high-entropy bubbles of expanding, neutrino-heated matter and 
are significantly more extended. The flux variations between 
low-emission and high-emission regions are usually 10-30%, 
but in very bright spots can reach up to about 70% for short 
times (typically of order milliseconds). These local variations 
are caused by enhanced fluxes of electron neutrinos and antineu- 
trinos, both being abundantly produced in the accretion layer, but 
the fractional variations are computed relative to the total neu- 
trino flux. The latter includes neutrinos of all three flavors, i.e. 
also heavy-lepton neutrinos and antineutrinos, whose emission 
(about 30-50% of the total flux) comes from the core and is not 
strongly afl'ected by accretion emission and thus more isotropic. 
After the explosion has been launched, typically a fewer accre- 
tion downflows survive and the number of hot spots decreases. 
The activity is now concentrated in a few regions (Fig.[9j middle 
panel) instead of showing the previous high-order multipole pat- 
tern that encompasses the whole sphere. After the accretion has 
finally ceased, the surface flux variations are determined by the 
convectively modified neutrino transport from below the neutri- 
nosphere. The granular pattern of the neutrino emission then 
reflects the smaller convective cells in the radially less extended 
convective layer inside the PNS (Fig. |9] lower panel). In this 
phase the variation amplitude of the total neutrino surface flux is 
usually around 5%, in peak activity periods up to 10%. 

During all three stages, however, the flux variation pattern is 
extremely time-dependent and non- stationary. The accretion hot 
spots appear, fade away again, and grow once more at difl'erent 
locations on the surface. Only in the phase of shock acceler- 
ation and simultaneous accretion, when the overall structure of 
the postshock ejecta begins to freeze out and a nearly self-similar 
expansion sets in, the coherence timescale of the emission pat- 
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Fig. 10. Left: Entropy distribution in a cross-sectional plane of model L15-1 at 1000 ms after bounce. The plane is chosen such that it contains 
the final vectors of the NS kick velocity induced by aspherical mass ejection (Vns, long white arrow pointing to the right) and of the kick velocity 
caused by anisotropic neutrino emission (Vns,v, short white arrow pointing towards the 7 o'clock position). The relative lengths of both arrows 
correspond to the absolute values of the velocities (Vns ~ 160km/s and Vns,v ~ 5 km/s), but the short arrow for Vns,v is stretched by a factor of 5 
for better visibility. A yardstick in the lower left corner of the figure indicates a spatial length scale of 2000 km. The contribution of anisotropic 
neutrino emission to the total NS kick velocity is clearly a small eff'ect. Right: Time evolution of the neutrino-induced kick velocity Vns,v and of 
the angle a^y between Vns and Vns,v. Note that the final velocities point roughly in opposite directions. 



tern increases. Nevertheless, there are still large temporal fluctu- 
ations because the accretion flows along the downflows are not 
continuous but change in strength. In the post- accretion phase 
the cellular pattern is uniform over the whole surface. The cell 
sizes show a tendency to increase at later times because the width 
of the convective layer grows with time. Simultaneously, the 
amplitude of the flux variations as well as the positions of cooler 
and hotter areas exhibit rapid temporal variability on timescales 
of milliseconds to some 10 milliseconds. 

As a consequence of the non- stationary, time-variable, and 
small-scale cellular structure of the surface-flux pattern, in par- 
ticular during the long-time neutrino emission in the Kelvin- 
Helmholtz cooling phase of the PNS when most of the neutrino 
energy is radiated away, we do not expect any sizable NS accel- 
eration due to anisotropic neutrino emission. 

Following [Scheck et al.| ( |2006| ), we calculate the neutrino 
momentum from the numerical results as 



Pv(0 ^ 



Jo Jr=Rob 



where py = IFyl/c"^ is the neutrino momentum density, and f 
is the unit vector in radial direction, which in our transport de- 
scription ("ray-by-ray approximation") also defines the direction 
of the neutrino flux (flux components in the angular directions 
are assumed not to exist). The first term accounts for the total 
momentum associated with all neutrino energy on the computa- 
tional grid at a given time t, while the second term is the time in- 
tegral of the neutrino momentum that has left the computational 
domain through the outer boundary until this time t. 



The numbers for the neutrino-kick velocity Vns,y listed for all 
models in Table [T] confirm our expectation: The data show that 
Vns,y makes only a very small, in most cases a negligible, correc- 
tion to the final NS kick velocity. It is typically much less than 
5 km s"^ except in model N20-l-lr, where Vns,y reaches 7 km s"\ 
which is still only ~4% of the mass-ejection induced kick veloc- 
ity Vns at 1.3 s after bounce. An entropy slice of model LI 5-1, 
which is a better resolved case with a relatively high neutrino- 
produced kick, is displayed at 1000 ms after core bounce in 
Fig. [T0| together with the time evolution of Vns,y and of the an- 
gle a^y between Vns,y(0 and Vns(0- The chosen slice contains the 
final (at ^ 1.3 s p.b.) vectors of both velocities. The final angle 
a^y is significantly larger than 100°, a result that is also found in 
most of our other 3D models (Table [T]). 

The velocity Vns,y increases steeply only after about 200 ms 
p.b. when a sizable accretion and emission asymmetry begins 
to develop. Before this phase the small angular variations of 
the neutrino emission do not allow for the growth of any no- 
^ Q-y^ ticeable value of Vns,y. Towards the onset of the explosion (at 

-400 ms; Table [T]) a nearly stationary pattern of accretion down- 
flows and high-entropy bubbles emerges in the postshock layer 
and persists during the subsequent expansion of the supernova 



shock and neutrino-heated ejecta. The steep rise of Vns,y between 
200 ms and 600 ms after bounce in model L15-1 is associated 
with a massive and extended, long-lasting downflow funnel that 
carries fresh gas swept up by the accelerating shock from the im- 
mediate postshock region towards the NS surface. This leads to 
a neutrino emission maximum in the right hemisphere (roughly 
at the 2:30 o'clock position) and thus to a neutrino-induced kick 
pointing in the opposite direction (short white arrow for Vns,y in 
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Fig. [To]). The relics of the massive downflow funnel are still 
visible to the right of the center in the left image of Fig. [10 



At about 600 ms p.b. the accretion ends and the neutrino emis- 
sion pattern adopts the high-multipole cellular structure char- 
acteristic of the PNS convection during the subsequent Kelvin- 
Helmholtz cooling phase (cf. above discussion and bottom panel 
in Fig. |9]). At this stage the neutrino-kick velocity has nearly 
reached its terminal value. It hardly changes later on because 
short-wavelength variations and short-timescale fluctuations di- 
minish the integrated neutrino emission asymmetry over longer 
periods of time. Nevertheless, a^y continues to exhibit some drift 
even between 600 ms and -1000 ms, because the gravitational 
pull by the anisotropic ejecta still causes an evolution of the as- 
sociated kick direction. Only at ^ > 1000 ms Vns has also found 
its final orientation about 135° away from Vns,y (Fig. [TO]), while 
the gravitational acceleration of the NS continues on a significant 
level even beyond 1.3 s after bounce (Table [T]). The neutrino kick 
typically damps the ejecta-induced kick: Large angles a^y are a 
natural consequence because long-lasting downflows determine 
the accelerating gravitational tug on the one hand and produce 
the dominant neutrino emission anisotropy on the other. While 
the NS kick associated with the asymmetric ejecta distribution is 
directed towards the downflow, the enhanced neutrino emission 
on the downflow side carries away momentum that leads to a NS 
recoil in the opposite direction. 

4. Neutron star spins 

Mechanisms that impart a natal kick to the NS can also lead to 
a NS spin even in the case of nonrotating progenitor stars. This 
is exp ected theoretically (Burrows et al.|1995t |Spmit & Phinney| 



1998[|Fryer & Young|2007 ) and suggested by observations (e.g. 
Farr et al.||2011| ). NS spin-up may, for example, be a stochas- 
tic process that is connected to random, nonradial impacts of 
the fast accretion downflows, by which low-entropy gas is chan- 
nelled through the gain region onto the PNS. Accretion down- 
flows, which hit the NS asymmetrically and not exactly head-on, 
exert torques and thus spin up the NS. A m ass Am co lliding with 
the NS with impact velocity Vi ~ Vff = ^/2GM^Jr and impact 
parameter d transfers an angular momentum of 



A/ns = Am 



2GMns 



6x 10 



46 



cm 

s 



Am_3 J30 Vi,io 



and for rigid rotation leads to a spin period of 



spin 



27r4 

A/ns 



0.2s(Am_3 J30 Vi,io) 



(18) 



(19) 



where Am is normalized by 10"^ M©, d by 30 km, vi is given in 
units of 10^^ cm s"\ and /ns ~ 2 x 10^^ g cm^ is the NS moment 
of inertia. 

Equations ([T8]) and ([T9|) give rough order of magnitude esti- 



mates and show that for moderate values of the nonaxisymmet- 
rically accreted mass (Am ~ 10"^ M©) the NS may receive an 
angular momentum of some 10^^gcm^s"\ corresponding to a 
specific angular momentum of order 10^^ cm^s"^ and spin peri- 
ods of some 100 ms. However, neither the number of accretion 
downflows nor their impact parameters can be computed analyt- 
ically. Quantitative answers require 3D hydrodynamical simula- 
tions of the postbounce accretion phase. 

Alternatively, it has been proposed that NS spin-up may not 
result from a stochastic process connected to random, nonradial 



impacts of accretion downflows. Instead, SASI spiral modes 
(m = 1,2 spherical harmonics components) might develop in 
the postshock accretion flow, by which a radial redistribution 
and separation of angular momentum between the compact rem- 
nant and the ejecta could be established: While NS spin is accu- 
mulated by the accretion of rotating gas, the ejecta carry away 
counter-rotating matter with an angular momentum of the same 
magnitude but with opposite direction (Blondin & Mezzacappa 



2007; Blo ndin & Shaw||2007t |Iwakami et al.| |2009; Fe rnandez 



|2010j . This phenomenon was seen in idealized numerical mod- 
els of collapsing stellar cores, making use of an ideal-gas EoS 
and ignoring neutrino heating (except by Iwakami et al. 2009]). 
It was also observed in laboratory experiments of a shallow water 
analogue of the SASI in two dimensions ( [Foglizzo et al.|2012 ). 
However, the timescale for SASI spiral modes to develop in the 
core without initial rotation was found to be close to one second, 
which is longer than the timescale for the onset of the explosion 
and for the saturation of the NS spins in all of our SN simula- 
tions. Correspondingly, we do not observe any clear signatures 
of SASI spiral modes in our models. In contrast, m = 1, 2 modes 
can dominate the asphericity of the shock surface, R^iO, cp), and 

of the postshock mass distribution, S(^, cp) = J^' drr^p (Eq. 3b, 
in models where the NS acquires only an unspectacular angular 
momentum. In all our models, however, the normalized am- 
plitudes of the m 9^ spherical harmonics components of the 
shock deformation never exceed a few percent (Fig. [5]). This is 
much smaller than the order unity shock displacements obtained 
by'Blondin & Mezzacappa' f2007)); [Blondin & S haw (2007) and 
^ernandez (2010). It is also possible that neutrino heating and 
neutrino-driven convection destroy the coherence of the spiral 
modes and thus impede their growth, or it might enhance mode 
coupling, thus amplifying energy redistribition between difl'er- 
ent modes. Moreover, even when shells with strong difl'eren- 
tial rotation in opposing directions arise in the turbulent region 
between PNS and shock, the question still remains how much 
of the induced angular momentum is finally accreted. [Rantsiouj 
|et al.| ( [2011) suspect that most of it is likely to be ejected in the 
explosion. 



4.1. Numerical results 

Because of violent convective activity, the matter in the neutrino- 
heated layer develops nonradial motion in diff'erent directions. 
These nonradial flows can be connected with sizable amounts 
of angular momentum. Integral values can indeed reach up to 
several 10"^^ gcm^ s"^ in difl'erent regions of the layer between 
PNS surface and accretion shock, compatible with the order-of- 
magnitude estimate of Eq. ([18). 

Since our simulations start from nonrotating collapsing stel- 
lar cores, any angular momentum acquired by the nascent NS 
must be balanced by angular momentum in the opposite direc- 
tion in the gas outside of the NS. Difl'erent from the linear mo- 
mentum of the compact remnant, which is calculated as the neg- 
ative value of the gas momentum integrated over the whole com- 
putational grid (see Eq.[T]), numerical tests showed that it is ad- 
visable to constrain the volume integration for the angular mo- 
mentum to the immediate vicinity of the forming NS. This helps 
minimizing numerical errors connected to the imperfect conser- 
vation of angular momentum of gas that expands outwards over 
radial distances extending over many orders of magnitude on a 
spherical grid with varied resolution. As in Wongwat hanarat] 



[etaLl ( [2QTOb 



we therefore determine an estimate of the angu- 
lar momentum Jns transferred to the NS by hydrodynamic forces 
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Fig. 11. Time evolution of the NS kick velocity (upper left panel) compared to the absolute value of the NS angular momentum (lower left 
panel), the polar and azimuthal angles (^spin, 0spin) of the NS spin direction (right top and middle panels), and the relative angle a^k between spin 
and kick directions (right bottom panel) for five selected models. Two of these models represent high-spin cases (N20-2, N20-3), two are low-spin 
cases (W15-6, W15-7), and one is an intermediate case (L15-3) from our model set. Note that the build-up of the angular momentum occurs in 
mainly a single impulsive episode during the convectively perturbed accretion phase of the NS. Together with its angles the angular momentum 
reaches its terminal value long before the NS acceleration has ended. Note that the clustering of all five models around a final ^spin ~ 90° is a mere 
coincidence, the whole set of computed models exhibits the expected random distribution in polar-angle space. 



associated with accreted or ejected gas from the following equa- 
tion: 



Jns(0 



WRns -'0 J47T 



where j is the specific angular momentum, Q denotes the solid 
angle, and Aq = Anr^ is the area of a sphere of chosen radius Vq. 
The first term in the brackets on the rhs of Eq. ( [2Q| ) accounts for 
the angular momentum contained by the spherical shell bounded 
by the NS "surface" at 7?ns (defined as the radius where the den- 
sity is 10^^ gcm"^) on the one side and bounded by Tq on the 
other. The second term gives the angular momentum that is car- 
ried by mass leaving the shell through surface Aq until time t. 
The introduction of the surface integral is possible in the con- 
text of NS spin-up, because angular momentum can be trans- 
ferred to the central object only when the latter exchanges mass 
with its environment. Long-range gravitational forces cannot ex- 
ert a torque on the compact remnant, because the NS (in a very 
good approximation) is described in our simulations as a gravi- 
tating point mass plus an essentially spherical near-surface layer 
of matter surrounding the inner grid boundary. The numerical 
evaluation of Eq. ( [2Q| ) yields results that do not depend on the 
exact location of Tq between 500 km and 1000 km. 

In TablefTlthe results for the angular momentum Jns, normal- 



the given numbers for Jns and the NS moment of inertia, /ns, 
as Tspin = 27r/ns/|JnsL Considering the NS as a rigid rotator and 
making the rough approximation that /ns is given by the Newto- 

(20) 

nian expression for a homogeneous sphere: /ns ~ |Mns^ns- "^^e 
(baryonic) NS mass Mns is also listed in Table [T] and we adopt 
a final NS radius of 7?ns = 12 km. Note that we assume that the 
angular momentum of the NS is conserved when the remnant 
contracts to its final radius during the neutrino-cooling evolu- 
tion. This ignores a possible -4 0% reduction associated with 
neutrino losses (see |Janka||2004| ) and with mass ejection in the 
neutrino-driven wind. 



ized by 10 g cm s" , the relative angle between spin and kick 
directions, c^sk, and the estimated final NS spin period, Tspm, are 
listed for all our model runs at the end of most of the simula- 
tions (1.1-1.4 seconds after bounce). We compute Tspm from 



The angular momentum values listed in Table [T] span a range 
from -0.1x10"^^ g cm^ s"^ to more than 7 x 10"^^ g cm^ s"\ while 
the corresponding spin periods vary by over a factor of 60, rang- 
ing from close to 0.1 s to nearly 8 s. Higher angular momenta 
usually lead to the smaller Tspm; this order is reversed by diflTer- 
ences of the NS masses only in a few cases. 

While there is no obvious systematic diflTerences of the NS 
angular momenta between the W15 and L15 models, although 
the L models clearly have the longest delays to explosion, three 
of the four members of the N15 series stick out having the high- 
est /ns values of the whole model set. In addition, two of the 
three B models are among those acquiring the lowest NS angu- 
lar momentum. These tendencies, cautiously interpreted in view 
of the limited number of models in each series, seem to be linked 
to the convectively stirred ejecta mass in the neutrino-heating re- 
gion. Assuming similar values for the specific angular momen- 
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turn j in all cases, the integration of Eq. (20) depends crucially on 
the mass involved. In particular model B 15 is special by its steep 
density decline outside the iron core in addition to the short evo- 
lution timescale to explosion, both of which tend to give a small 
angular momentum transfer to the accreting remnant. 

The values of the final angle between spin and kick vectors, 
Qfsk, in Table [T] do not reveal any clear correlations beween spin 
and kick directions, in particular no strong tendency towards 
alignment or counter-alignment. The distribution is nearly sym- 
metric for a^sk < 90° and a^sk > 90° (9:11 cases) and essentially 
flat when sorted into 60° intervals (7:6:7 cases), although a sta- 
tistically equal distribution in all directions would require twice 
as many events in the equatorial belt than in each polar wedge. 
Our model sample, however, is still too small to draw firm con- 
clusions in a statistical sense. 



4.2. Neutron star spin-up mechanism 

It is important to note that kicks and spins are attained by the 
NS through diff'erent mechanisms in our simulations, although 
both phenomena are connected to the development of large-scale 
nonradial flows in the collapsing stellar core, which lead to mo- 
mentum and angular momentum exchange between the accreting 
compact remnant and the asymmetrically expanding ejecta gas. 
In order to understand the difl'erences, we present in Fig. [TT] the 
time evolution of /ns compared to Vns for a selection of represen- 
tative models, namely for the two high-resolution N15-runs with 
the largest values of /ns, for the two W15-cases with the lowest 
/ns, and for model LI 5-3 as an intermediate case. The left two 
panels show that the angular momentum increases much faster 
and reaches a value close to its asymptotic one much earlier than 
the recoil velocity does. Difl'erent from the NS kick most of the 
NS spin is attained in an impulsive event within a time interval 
of 1 00-200 ms relatively early after core bounce. No significant 
/ns evolution takes place beyond t > 500 ms, i.e. after the onset 
of the explosion (measured by ^exp). In contrast, the NS velocity 
is very small at the time when the PNS receives its angular mo- 
mentum, but the kick can grow monotonically with considerable 
rate for several seconds. 

The reason for this diff'erence is the fact that angular mo- 
mentum can be transferred to the NS only by nonspherical mass 
flows and thus during a few hundred milliseconds of accretion. 
This takes place before the onset of the explosion and can go on 
afterwards only as long as accretion still continues simultane- 
ously with the outward expansion of the SN shoc]^ In contrast, 
the recoil of the newly formed NS is mainly established by the 
long-range anisotropic gravitational attraction between the com- 
pact remnant and the asymmetric ejecta as described in detail in 
Sect. |3] In Fig. 1 1 this difl'erence manifests itself in the long- 
lasting, monotonic rise of Vns(0 during the whole displayed time 
and even much longer (cf. Fig. [4]), whereas most of the angular 
momentum is accumulated in the NS during a period of mas- 
sive accretion shortly after the convective and SASI instabilities 
have grown to their full strength. When the accretion to the NS 
is diminished by the outward acceleration of the gas behind the 
explosion shock, sporadic downdrafts reaching the NS can still 
lead to modest variations of /„« for some 100 ms. The duration 



Hi mass per 
solid angls [g] 




Note that the neutrino-driven wind that develops after the onset of the 
explosion is essentially spherically symmetric in our simulations. Pos- 
sible angular variations caused by PNS convection have low amplitudes 
and short wavelengths and thus do not carry any significant amount of 
momentum. Also the rate of angular momentum loss in the dilute out- 
flow is typically very low. 



Fig. 12. ^^Ni distribution for model W15-2 at 365 ms and 515 ms 
after bounce. The integrated nickel mass per solid angle is plotted in 
the cross-sectional plane that contains the NS kick (indicated by the 
black arrows pointing to the right) and spin vectors. Red colored regions 
show that more nickel is formed in the direction opposite to the NS kick 
direction. The size of the displayed volume can be estimated from the 
yardsticks given in the lower left corner of each image. 



of this phase depends on the progenitor structure and the detailed 
dynamics of the explosion, but even for relatively low explosion 
energies it does not last any longer than -500 ms after the revival 
of the SN shock. 

The difl'erences between the kick and spin mechanisms are 
reflected not only by the evolution of the NS kick and spin mag- 
nitudes but also by the evolution of the kick and spin directions. 
During the initial steep growth of /ns the spin direction (^spin, 
0spin) mostly settles to its final direction, in particular in the two 
high-spin N20-models. Long-time variations of the NS angular 
momentum then correlate with only smaller changes of the spin 
direction, which signal the impact of new accretion downdrafts. 
Such changes become less frequent at later times and finally die 
away when the accretion ends. In the L15-models, which have 
the highest densities in the shells surrounding the iron core and 
therefore explode with the longest delays (a good example is 
model L15-3, where the shock is revived as late as -500 ms af- 
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ter bounce), the period of large /ns variations is more extended. 
Slow fluctuations of the spin- vector angles and a slight drift of 
/ns occur even until more than one second post bounce because 
of the angular momentum transported away with the more mas- 
sive neutrino-driven wind in the L models, in contrast to the W 
and N models. The long-time variation of the spin difl'ers from 
that of the kick, whose direction freezes out as soon as the bulk 
of the mass carrying the ejecta asymmetries expands essentially 
self- similarly, although Vns continues to increase. The bottom 
right panel of Fig. [TT] provides the time evolution of the relative 
angle between kick and spin vectors. The high-frequency vari- 
ations of Qfsk are caused by rapid changes of the kick direction 
(^kick, 0kick) until about 200 ms after the onset of the explosion, 
while a slow drift of a^sk on longer timescales (up to ~ 1 s in the 
most extreme case) is connected to an ongoing reorientation of 
the spin direction. 

Interestingly, the highest NS kicks and spins are not ob- 
tained for the same models or progenitor. While the remnants 
of the N20 progenitor exhibit the tendency to develop the high- 
est spins, those of the W15 series are found to end up with the 
largest kicks. This can be considered as another manifestation 
that spin and kick mechanisms are linked to difl'erent aspects 
of anisotropic mass flows in the SN core, namely the former to 
asymmetries of the accretion flow and the latter mainly to as- 
phericities of the ejecta distribution. 

It is unclear whether NS spin periods that are more extreme 
than those found in some of our N20 models, in particular pe- 
riods below 100 ms, can be obtained in nonrotating stars. It 
is possible that faster NS rotation requires the progenitor core 
to rotate (e.g., Heger et al.| |2005[ |Qtt et al.| |2006). In a rotat- 
ing environment, however, the growth of SASI spiral modes will 
be severely altered and potentially fostered ([Blondin & Mezza-| 



cappa||2007[|Yamasaki & Foglizzo [20081 [Iwakami et al.||2009| ), 
and three-dimensional simulations are therefore indispensable to 
make predictions for the connection between progenitor and NS 
rotation. 

5. Heavy-element production and explosion 
asymmetries 



Fryer & Kusenko| ^006*) pointed out that ejecta asymmetries 
could be used to distinguish between difl'erent theoretical sug- 
gestions for the kick mechanism. Ejecta-driven kicks like the 
gravitational tug-boat mechanism discussed here must be ex- 
pected to be associated with stronger ejecta motion, i.e. with 
higher ejecta momentum, in the direction opposite to the kick. 
However, diagnosing momentum asymmetries of the ejecta may 
be difficult except in cases where the center of mass of the 
gaseous SN remnant is clearly displaced relative to the NS posi- 
tion or where the gas cloud exhibits extreme deformation with- 
out gradients of the environmental conditions being responsible 
for the asymmetric ejecta expansion. The SN remnants Cas- 
siopei a A ([Isensee et al. 2010; DeLaney et al. 2010; Rest et al.| 
20TTt [Hwang & Laming 2012J andPuppis A (^Petre et al.|1996[ 



Winkler & Petre|2007t [Katsuda et al.|2008||2010t [Becker et al.| 
20T2V may be such lucky cases, because in the former all the 
SN ejecta in the remnant, whose mass is dominated by oxygen, 
while in the latter fast-moving, oxygen-rich optical filaments and 
knots seem to have been recoiled opposite to the direction of the 
high-velocity compact stellar remnant. In general, however, a 
major problem for observing momenta asymmetries linked to 
an asymmetric beginning of the explosion may result from the 
fact that the outgoing SN shock transfers the bulk of the explo- 
sion energy and ejecta momentum to the He and H shells, which 



dominate the ejecta mass in events other than stripped SNe. An 
initially asymmetric explosion thus becomes much more spher- 
ical when the SN shock sweeps up massive stellar He and H 
layers. In this context it is important to note that the NS kinetic 
energy and the NS momentum are very small compared to the ki- 
netic energy of the SN explosion, £^k,SN. and the corresponding 
measurable momentum, Pk,SN = y2Mej£'k,sN, associated with 
the radial motion of the ejecta mass Mej. To determine the (rela- 
tively small) linear momentum of the ejecta, which balances the 
momentum of the NS, it is necessary to sum up all components 
of the momentum vectors of the whole ejecta mass. Because of 
projection effects this may be difl&cult and error-prone even if all 
of the ejecta can be observed. 

For this reason we propose here to look for correlations be- 
tween NS kicks and asymmetries in the distribution of heavy ele- 
ments. These are explosively produced in the innermost, shock- 
heated (and neutrino-heated) SN matter during the first seconds 
of the blast, and their globally asymmetric distribution will not 
be destroyed lateron. In particular iron-group nuclei, most of 
which are radioactive ^^Ni decaying to stable ^^Fe, are good trac- 
ers of such early explosion asymmetries. This is connected to the 
fact that the initial Fe-group material in the progenitor core is 
completely photo-disintegrated during collapse and then buried 
in the forming NS. All expelled Fe-group nuclei are freshly as- 
sembled in the neutrino-heated SN ejecta and in the silicon and 
oxygen layers that are heated to sufficiently high temperatures 
by the accelerating SN shock. 

Figur e visualizes the corresponding directional asymme- 
try of the ^^Ni production in the plane containing the NS veloc- 
ity and spin vectors in the high-kick model W15-2. The images 
show the nickel mass per unit solid angle at post-bounce times 
of 365 ms and 515 ms (the explosion in this model sets in about 
250 ms after bounce; Table [T]), which is obtained by integrating 
the nickel-mass density over radius within each zone: 



rri+m 



drrV(r)XM(r), 



(21) 



where XNi(r) is the mass fraction of nickel and p(r) the mass 
density at location r. The color coding clearly shows more nickel 
production (intense red) in the direction opposite to the kick vec- 
tor. 

At around 450-500 ms after bounce the nickel fusion in the 
shock-heated layers comes to an end. While most of the nickel 
expands in a thin and dense shell behind the shock, some of 
it is channeled into massive clump-like "pockets" between the 
more dilute high-entropy bubbles of rising and outward pushing 
neutrino-heated matter. Also these bubbles may contain a fair 
amount of nickel from the a-rich freeze-out that occurs during 
the expansion cooling of the bubbles. However, because these 
regions have a slight neutron excess in our models, the nuclear 
flow is assumed to lead to the production of the "tracer" species 
included in our network as described in Sect. 12.31 Since the 
actual composition of this tracer material strongly depends on 
the inaccurately determined value of set by and absorp- 
tion and emission in the neutrino-processed ejecta, the amount of 
^^Ni in this material is uncertain and we have not depicted it in 
Fig. 12 The tracer material, however, enhances the hemispheric 
asymmetry of the nickel distribution. This is visible from the 
data listed in Table [2j which provides abundance yields of the 
nuclear species in our burning network after all nucleosynthesis 
processes have ended (depending on the progenitor and explo- 
sion energy, this is the case typically 100-300 s after the onset 
of the explosion). Higher values of ^^Ni in the southern hemi- 
sphere correlate with a similarly large excess of tracer material 
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Table 2. Hemispheric ejecta yields of helium and metals, and total nickel mass for the high-kick models W15-1 and W15-2, and the moderate- 
kick models L15-1 and L15-2. The numbers are integral values of the elemental abundances in the whole stellar matter that gets expelled by the 
supernova explosion. The tracer column provides the yield of Fe-group nuclei in neutrino-processed ejected material, some undetermined fraction 
of which may be ^^Ni. The north-polar direction is defined to coincide with the NS kick direction. 



Model 


'He [Me] 


[10- 


'Mo] 


[10- 


■'Mo] 


20Ne [10 


-'Mo] 


24Mg 


[10-2 Mo] 




North 


South 


North 


South 


North 


South 


North 


South 


North 


South 




W15-1 


2.78 


2.66 


1.18 


1.10 


3.68 


3.75 


8.90 


8.49 


2.41 


2.85 




W15-2 


2.78 


2.65 


1.16 


1.12 


3.43 


3.84 


8.67 


8.49 


2.16 


2.86 




L15-1 


2.39 


2.34 


0.90 


0.87 


2.77 


2.89 


5.00 


5.06 


2.12 


2.49 




L15-2 


2.40 


2.39 


0.89 


0.87 


2.85 


2.79 


5.21 


4.88 


2.47 


2.42 






Model 


28Si [10- 
North 


-2 Mo] 
South 


40Ca [10 
North 


-2 Mo] 
South 


44j- [10- 
North 


"'Mo] 
South 


^^Ni [10- 
North 


-2 Mo] 
South 


Tracer [10"^ Mo] 
North South 


Total ^^Ni Mass [Mo] 


W15-1 


1.88 


2.92 


1.33 


4.81 


0.68 


2.43 


1.26 


4.28 


2.23 


6.08 


0.055-0.139 


W15-2 


1.74 


2.83 


1.27 


4.66 


0.81 


2.17 


1.37 


4.09 


2.22 


6.27 


0.055-0.139 


L15-1 


1.75 


2.33 


1.76 


2.47 


1.49 


2.40 


1.34 


1.87 


4.78 


7.20 


0.032-0.152 


L15-2 


2.13 


2.15 


2.54 


2.74 


2.32 


2.55 


1.81 


1.89 


8.68 


9.74 


0.037-0.221 



on the same side. We define the northern hemisphere as the half- 
sphere around the NS kick direction, which implies that more 
^^Ni and tracer material are ejected opposite to the NS kick di- 
rection. This eflTect is especially pronounced (with hemispheric 
contrast up to a factor 3.5) in high-kick models like W15-1 and 
W15-2. The last column of Table |2] provides our determination 
of the total ^^Ni yields of the listed models. The uncertainty 
range is defined by the assumption that in the lower limit the 
tracer material contains no ^^Ni and in the upper limit all tracer 
material is ^^Ni. 

The comparison of the high-kick models W15-1 and W15- 
2, both of which have NS recoil velocities in excess of 
500 km (Table[l]), with the moderate-kick model LI 5-1 (Vns ~ 
250km s'l) and model LI 5-2, whose NS receives a recoil ve- 
locity of only ~100kms-^ is enlightening. The data in Table [2] 
demonstrate that the hemispheric diflTerences of nickel and tracer 
ejection are largest for W 15-1 and W15-2, still considerable for 
LI 5-1, and nearly disappear for LI 5-2. Similarly large hemi- 
spheric diflferences as for ^^Ni and our tracer species show up in 
the masses channeled by the network into ^^Si, ^^Ca, and ^^Ti. 
All of these nuclear specie^ emerge from the innermost ejecta 
regions, which carry the explosion asymmetries and the mass- 
distribution inhomogeneities that are responsible for the gravita- 
tional acceleration of the newly formed NS. While some degree 
of hemispheric diflTerences correlating with the NS kick magni- 
tude is still visible for ^^Mg, the eflfect in the case of ^^O is very 
small. Although oxygen is more abundantly synthesized behind 
the stronger shock in the southern hemisphere, an asymmetry of 
the distribution of this element is largely obscured by the dom- 
inance of the spherically symmetric ejecta mass that originates 
from the unperturbed oxygen shell of the progenitor star. The 
other lighter species (^^Ne, ^^C, and "^He) do not exhibit the 
southern excess. In contrast, their yields are slightly reduced 
in the southern hemisphere — most pronouncedly in the case of 



^ It should be noted that the small number of species in the employed 
network is likely to be responsible for a considerable overestimation of 
the "^"^Ti production in our simulations. A better determination of the 
^Ti yields will require more sophisticated network calculations. 



2^Ne — because of the explosive production of the higher-mass 
nuclei at the expense of these lighter abundances. 

The dependence of the explosion strength on direction is vi- 
sualized in Fig. [13] in terms of the variations of the shock radius 
and thus shock velocity at around the time when the nickel nu- 
cleosynthesis has just taken place. Figure [14] shows the corre- 
sponding radially integrated nickel masses per solid angle. 



dMi 



Ni 



dQ 



(0, 



drrV(r)XNi(r), 



(22) 



where the radial integration is performed between the inner and 
outer grid boundaries, 7?ib and 7?ob, respectively. In Fig.pT 



we 



provide volumetric impressions of the nickel distribution in the 
ejecta by plotting isosurfaces of a suitably chosen constant value 
of the nickel mass per grid cel|^ The coordinate system in these 
images has been rotated such that the NS kick vector points to 
the 12 o'clock position. 



In Fig. 13 one can see clear diflTerences of the shock defor- 
mation between low-kick and high-kick models. In the high-kick 
cases of W15-1 and W15-2, the northern side, which is defined 
as the hemisphere around the NS kick direction (marked by the 
red cross), contains a huge, deep trough (blue and black) whose 
minimum roughly coincides with the kick direction. In contrast, 
the southern hemisphere is filled by a wide bump with several lo- 
cal maxima. Consequently, there is a clear dipolar shock defor- 
mation along the kick axis (cf. also Fig. [5]). The moderate-kick 
model LI 5-1 still exhibits a tendency for the same asymmetry 
pattern with a shock minimum close to the red cross and a max- 
imum near the white plus sign in the opposite direction. Now, 
however, shallow minima and maxima also appear in the south- 
ern and northern hemisphere, respectively, indicating the growth 
of a quadrupolar deformation component. The low-kick model 
LI 5-2 violates the north-south imbalance. Both shock maxima 



^ The displayed quantity has the advantage that the bulk of the pro- 
duced nickel becomes visible. Both the nickel density, which essentially 
follows the radial gradient of the mass density, and the nickel mass frac- 
tion, which does not necessarily trace regions with high integral nickel 
mass, are less suitable for this purpose. 
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W15-1 




Fig. 13. Shock deformation for the high-kick models W15-1 and W15-2 and for the moderate-kick models L15-1 and L15-2 at 1.3 s or 1.4 s 
after bounce, as given in the upper right corner of each map. These times are close to the evolution phases when the nickel production in the 
shock-heated ejecta takes place. The color scale is linear with black indicating minima and red maxima. The differences between maxima and 
minima relative to the average shock radii are 30%, 35%, 34%, and 29%, respectively. Red crosses indicate the NS kick directions while white 
plus signs mark the directions opposite to the kicks. The white lines separate the corresponding two hemispheres. The northern hemisphere is 
defined as centered around the NS kick vector. 



and minima cross the equatorial line and are obviously not cor- 
related with the locations of the cross and plus signs. Instead 
of the dipolar shock deformation of models W15-1 and W15-2 
we now see a quadrupolar pattern, consistent with the left lower 
panel of Fig. [5] We note that the relative variation of the shock 
radius in all four displayed cases is similar. The shock deforma- 
tion in high-kick and low-kick models therefore does not need to 
differ in the variation amplitude, but instead the type of dominant 
spherical harmonics mode can be the more important factor. 



In Fig. 14 the highest concentrations of ejected nickel are 



found exactly at the locations of the shock maxima in Fig. [TT 
while big "holes" of the nickel distribution appear in the regions 
of shock minima. Again, the angular pattern of the ^^Ni distribu- 
tion of the low-kick model LI 5-2 is clearly dominated by higher 
spherical harmonics modes than in the high-kick models W15-1 
and W15-2. Model LI 5-1 with a moderate kick defines an inter- 
mediate case. The directional asymmetries depicted in Fig. [14 
are complemented by the volumetric information of Fig. [15 The 
black areas seen in Fig. [14] correspond to directions where holes 
exist in the spatial nickel distribution, while the green-yellow- 
red maxima in Fig. [14] appear as big, clump-like, coherent struc- 
tures of nickel in Fig. 15 The hemispheric nickel imbalance 
in models W15-1 and W15-2 is obvious, whereas L15-1 shows 



a bipolar structure that is misaligned with the NS kick vector. 
Dififerent from the three others, LI 5-2 exhibits a more isotropic 
nickel distribution and considerably more power in structures on 
smaller angular scales. 

In summary. Figs. [T2}|T5] confirm that the production of 
nickel and of adjacent products of explosive burning (cf. Table[2]) 
are excellent tracers of explosion asymmetries. An observational 
determination of major anisotropics of the motion of heavy ele- 
ments between silicon and nickel in SN remnants and an expan- 
sion of these metals mainly opposite to the direction of the NS 
motion would lend strong support to the NS acceleration by the 
described gravitational tug-boat mechanism. 



6. Summary and conclusions 

We have performed a set of 20 SN simulations from shortly after 
core bounce to more than one second after the onset of the explo- 
sion for four diflTerent 15 M© and 20 M© nonrotating progenitor 
stars, making use of an axis-free Yin- Yang grid in three dimen- 
sions. The grid helped us to avoid numerical artifacts as well as 
severe time- step constraints near the polar axis of a usual spheri- 
cal grid. A further gain in performance was achieved by excising 
the high-density core of the PNS (at densities significantly above 
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Fig. 14. Integrated nickel masses per solid angle in all directions displayed in 4?: maps for models W15-1, W15-2, L15-1, and L15-2 at 
postbounce times well after all nucleosynthesis processes have seized in our simulations. Red crosses indicate the NS kick direction while white 
plus signs mark the direction opposite to the kick. The corresponding hemispheres are separated by a white equator line. Note the differences 
between the high-kick W15 models and the moderate-kick LI 5 models. The former show a clear enhancement of the nickel production in the 
hemisphere opposite to the kick, whereas the latter eject nickel significantly more isotropically. 



the neutrinosphere) and replacing it by an inner boundary condi- 
tion, where the contraction of the PNS core was imposed and the 
neutrino luminosities were chosen such that neutrino-energy de- 
position behind the shock led to SN explosions of defined energy. 
The neutrino transport on the grid was treated by a computation- 
ally efiftcient, analytic characteristics integrator, in which a grey 
approximation by was used averaging over Fermi-Dirac spectra 
and flux variations with latitude and longitude were represented 
by employing a ray-by -ray description (see |Scheck et al.|20Q6| ). 

A subset of our explosion models was evolved beyond three 
seconds postbounce to follow the approach of the NS recoil ve- 
locity to the saturation level. Some of these simulations were 
carried on for several 100 seconds to determine the explosive 
nucleosynthesis of heavy elements and the spatial distribution of 
the burning products. 

Our results demonstrate the viability of the gravitational tug- 
boat mechanism for NS acceleration in three dimensions. In this 
scenario the newly formed NS is accelerated over a timescale 
of many seconds mainly by the anisotropic gravitational attrac- 
tion of the asymmetrically expanding neutrino-heated ejecta gas, 
in which hydrodynamic instabilities had created large-scale as- 
phericities prior to the onset of the explosion. Before the ac- 
cretion flows to the NS are quenched by the accelerating explo- 
sion, hydrodynamic forces on the NS play a nonnegligible role. 



too. We emphasize, however, that anisotropics of the neutrino 
emission, which our ray-by-ray transport approximation tends 
to overestimate, contribute to the NS kick velocity only on a mi- 
nor level of percents or less. The reason for this small eflTect 
is simple: Anisotropics of the neutrino emission exhibit short- 
timescale intermittency and are characterized by a higher-order 
multipole structure in angular space. As a consequence, the 
neutrino-induced momentum transfer is diminished by statisti- 
cal averaging and is therefore unable to cause any significant net 
kick in a certain direction. Moreover, only a small fraction of the 
total neutrino energy loss is radiated anisotropically during the 
simulation time. Instead, the far majority of the escaping neu- 
trinos contribute to the isotropic background flow that diflTuses 
out of the essentially spherical, hot accretion layer around the 
nascent NS. 



NS acceleration by the gravitational tug-boat mechanism 
was first discussed on the basis of axisymmetric (2D) simula- 
tions by |Scheck et al.| ( |2004[ |2006| ) and was confirmed in 3D 
with a small set of explosion models by |Wongwathanarat et al. 
(2010b). Moreover, support for this mechanism was recently 
also provided by the 2D simulations of Nordhaus et al. (2010, 
2012| ), in which the NS was allowed to move self-consistently 
out of the grid center. 
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Fig. 15. Volumetric three-dimensional visualization of the nickel distribution for models W15-1, W15-2, LI 5-1, and LI 5-2 at the postbounce 
time given in the top left corner of each panel along with the model name. The semi-transparent isosurfaces correspond to a chosen value of the 
nickel mass per grid cell of 3 x 10^^ g and are displayed at a stage well after all nucleosynthesis processes have seized in our simulations. The 
orange vectors represent the NS kick directions and are scaled by the corresponding NS velocities (96 km s~^ for the shortest arrow and 575 km s~^ 
for the longest one). The high-kick models W15-1 and W15-2 in the upper two panels exhibit a clear asymmetry with much more nickel being 
ejected in the hemisphere opposite to the kick direction. In contrast, the moderate-kick models LI 5-1 and LI 5-2 in the lower two panels exhibit 
a more isotropic distribution of the nickel, in particular no obvious hemispheric asymmetry between kick and anti-kick directions. Note that the 
radial distribution of the nickel may be strongly affected and changed by subsequent mixing instabilities that develop after the outgoing shock has 
passed the composition-shell interfaces of the progenitor star. The hemispheric differences in the nickel ejection, however, will not be destroyed 
during the later supernova explosion. 



While the artificial constraint to axisymmetry enforces a col- 
limation of large-scale flows parallel to the polar grid axis and 
thus tends to favor a pronounced, dipolar deformation of the ex- 
plosion, the 3D asphericities appear less extreme and seemingly 
less promising for high NS kicks. Nevertheless, the inhomo- 
geneities and anisotropics in the massive postshock shell are siz- 
able also in 3D. In fact, low spherical harmonics modes of dipo- 
lar and quadrupolar character can dominate the asymmetry of 
the shock and postshock ejecta also in 3D, although with some- 



what smaller amplitude than in extreme 2D situations (Fig. [5]). 
Correspondingly, also in our set of 20 three-dimensional mod- 
els several cases were found that produced NS kick velocities 
well above of 500kms"\ and our current record holder devel- 
oped a NS recoil of more than 700 km s"^ after 3.3 s (Table \\\ 
and Fig. [^. This is compatible with the 2D results of Scheck] 



et al. 



(2006), who obtained velocities up to 800 km s ^ after one 



second and estimated the final NS speed in one out of more than 
70 models to be nearly 1000 km s"^ and in another one to be 
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more than 1200 km s"^ Simple analytic estimates of the long- 
time influence of the gravitational forces by the expanding SN 
ejecta confirm that dense clumps of some 10"^ M© (accounting 
for mass distribution asymmetries of around 10% in the post- 
shock ejecta shell) could well accelerate the NS to velocities 
of 1000-2000 km s"^ as measured for the fastest young pulsars 
(cf. Eq. 10). With an ejecta-shell mass of a few 10"^ M©, a 



similar accelerating eff'ect on the NS can be achieved by hemi- 
spheric asymmetries of the expansion velocity of several 10% 
(cf. Eq. [12]), which is in the ballpark of shock expansion diff'er- 
ences seen in some of our models (Fig.fTS ). 

The long-time acceleration of the NS is therefore connected 
to slow-moving, massive ejecta "clumps", whose expansion lags 
behind the rest of the ejected material of the SN core. For 
this reason the NS gains momentum opposite to the direction 
of fastest expansion of the explosion shock and postshock shell. 
Correspondingly, the explosive nucleosynthesis of heavy ele- 
ments from ^^Si to the iron group, especially also of radioac- 
tive ^^Ni, occurs preferentially in the hemisphere pointing away 
from the NS velocity vector. This eff'ect clearly correlates in 
strength with the size of the kick. We therefore propose that 
an observational confirmation of opposite momentum directions 
of the NS and of the bulk of elements heavier than silicon in 
supernova remnants would lend support to the gravitational kick 
mechanism discussed in this paper. Such an analysis for the Cas- 
siopeia A (Cas A) SN remnant with its central point source could 
be enlightening. 

Nonradial (off'-center) accretion downflows in the postshock 
layer can also lead to a spin-up of the nascent NS in nonrotating 
progenitor stars. This requires angular momentum separation 
in the postshock flow during the postbounce accretion phase of 
the PNS. The corresponding angular momentum can accumulate 
to values of several 10^^ gcm^ s"^ in our simulations, in agree- 
ment with an order-of-magnitude estimate on grounds of simple 
considerations (Eq.fTSl). Assuming angular momentum conser- 
vation during the later evolution, we derive spin periods between 
~100 ms and several seconds for the NS in our models. We note, 
however, that angular momentum conservation is difficult to as- 
sure in the numerical models, and thus these numbers should be 
considered only as a rough range of possibilities. 

While the NS acceleration is a long-time (continuing over 
seconds) phenomenon mostly achieved by the gravitational pull 
of the postshock ejecta shell, the spin-up of the nascent NS is 
connected to anisotropic mass flows that transfer angular mo- 
mentum to the accretor. The relevant postbounce accretion phase 
of the PNS may last for several 100 ms beyond the onset of the 
SN blast. This limits the time interval for the spin-up of the NS 
to a much shorter period than the NS acceleration. Kick and spin 
of the NS in this scenario are therefore not obtained by the same 
process: The former is linked to mass ejection asymmetries, the 
latter to anisotropic, nonradial accretion flows; they are mediated 
by difl'erent forces — gravitational versus hydrodynamical — , 
and they also build up on difl'erent timescales. It is therefore 
not astonishing that we do not find any correlation of spin and 
kick directions in our explosion models of nonrotating progeni- 
tor stars, and a possible spin-kick alignment suggested by obser- 
vat ions (e.g., [Jo hnsto n et al. 2 005 , 2007 ; Kaplan et al. 2008 ; also 
see [Wang et a l. 2006 , 2007] ) remains unexplained. 

In this context, the question of the consequences of (even 
a modest amount of) angular momentum in the progenitor core 
is interesting and deserves future exploration. In the described 
kick scenario, however, rotational averaging of the kick, which 
could cause an alignment of pulsar velocity and spin as sug- 
gested by |Lai et al.] ( |2001| ) and |Wang et aL] ( |2007| ), is unlikely 



to be efl&cient despite the long timescale of the NS acceleration. 
Rotational motions in the expelled gas, whose clumpy asymme- 
tries accelerate the NS by gravitational forces, are diminished 
during expansion, because under the constraint of angular mo- 
mentum conservation the angular velocity decays with growing 
radius like r"^. It is therefore not possible to destroy the ran- 
domness of the kick direction just by the rotation of the accel- 
erating ejecta asymmetries. An asphericity that preferentially 
accelerates the NS along the spin direction would have to be im- 
printed in the ejecta gas already at the onset of the explosion. It 
is unclear whether core rotation can produce such an asymme- 
try in the case of reasonable and realistic assumptions about the 
spin rate of the progenitor core prior to colla pse. With predicted 
typical Fe-core rotation periods of (9(100 s) ( Heger et al.|2005 ), 
which lead to NS spin periods of (9(10 ms), the neutrino-heated 
postshock layer has rotation periods of about one second. This 
is much too slow to have any dynamical influence. On the other 
hand, if sufl&ciently fast rotation allowed for a rapid growth of 
SASI spiral modes or triaxial dynamical or secular instabilities, 
kicks perpendicular to the NS spin appear more likely than the 
alignment of both. In our simulations we do not find any evi- 
dence for the large-amplitude SASI spiral modes that were ob- 
served in idealized setups (without neutrino eff'ects) by ( ,Blondin] 
|& Mezzacappa|20()7]) an d [Fernandez| ( |2010| ) and experimentally 
by |Foglizzo et al.H2012| ). The growth of such large-scale, coher- 
ent patterns in the postshock accretion flow might be inhibited by 
the presence of neutrino heating, or their growth timescale might 
be longer than the evolution period until the onset of the explo- 



sion in our models. Numerical studies ( Blondin & Mezzacappa 
2007 [ llwakami et ^^1200 9^) and analytical treatment (Y amasaki 



& Foglizzo 2008 ), however, predict faster growth of SASI spiral 



modes in rotating environments. 

Our still limited set of simulations does not allow us to draw 
any firm conclusions on progenitor or explosion dependences of 
the spin and kick magnitudes. One of the employed 15 M© stars 
(the W15 progenitor) seems to exhibit a tendency to favor higher 
recoil velocities on average than the other three investigated pro- 
genitors, but the exact reason for this trend is not determined. 
Podsiadlowski et aL] ([2004]) hypothesized that the collapse 



and explosion of low-mass SN progenitors with ONeMg cores 
leads to NS with low- velocity kicks. They argue that in the case 
of ONeMg cores in contrast to more massive Fe-core progen- 
itors the explosion sets in more quickly and does not provide 
the time for the growth of large- amplitude asymmetries of low- 
order multipole type in the neutrino-heated accretion layer be- 
tween NS and stalled shock. Indeed, the bounce shock acceler- 
ates quickly in the steep density profile at the edge of collaps- 
ing ONeMg cores and the explosion sets in on a timescale of 
only- 100 ms after core bounce (Kitaura et al.|2006]|Janka etaL] 



2008 ). This leaves time for only one convective turnover. There- 
fore, the angular size of the convective plumes remains rather 
small. Because of the extremely fast, spherical expansion of the 
shock, the convective structures are also rather symmetric in all 
directions (cf. Fig. 1 in |Wanajo et al.|[2011| ). These facts to- 
gether with the extremely small mass swept up from the dilute 
environment of the ONeMg core by the outgoing shock disfa- 
vor high NS kicks by the gravitational tug-boat mechanisnj^ In 
our scenario of NS acceleration the diff'erences of low-kick and 
high-kick cases are connected to diff'erences in the explosion dy- 
namics, manifesting themselves in a diff'erent growth of hydro- 



We note in passing that 



Podsiadlowski et al. (12004 as well as 



|den Heuvel ( 2010 ) erroneously attributed the acceleration of the NS to 
anisotropic neutrino emission instead of asymmetric mass ejection. 
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dynamic (convective and SASI) instabilities during the onset of 
the explosion. The growth of these instabilities can be seeded 
by small random perturbations (which should be present in any 
convectively perturbed stellar environment and are imposed by 
us artificially on spherical pre- or post-collapse progenitor mod- 
els) and leads to anisotropic mass ejection. Large-scale asymme- 
tries that might develop in the progenitor core during convective 
shell-burning phases prior to collapse as advocated by Burrows 



& Hayes| ( [T995| ); |van den Heuvel| ( [2QT0l ), and |Arnett & Meakin 
^ 2Q11| ) are therefore an interesting possibility but do not seem to 



be necessary for an explanation. Of course, if such large-scale 
asymmetries of the progenitor structure occurred, they could 
have important consequences for the SN explosion and the NS 
kick mechanism. Unquestionably, the availability of more realis- 
tic, three-dimensional progenitor models for core-collapse sim- 
ulations would be highly desirable. 
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